GAMs and GPs

import pandas as pd
import matplotlib.pyplot as plt
import bambi as bmb
import seaborn as sns
from pygam.datasets import mcycle
from pygam import LinearGAM, s, f, GAM, l, utils
import numpy as np
import arviz as az
import pymc as pm


random_seed = 100

import warnings

warnings.simplefilter('ignore')

Flexibility and Calibration

In this blog post we’re going to dive into modelling of non-linear functions and explore some of the tooling available in the python eco-system. We’ll start by looking into Generalised Additive Models with splines in pyGAM before preceding to look Bayesian versions of spline modelling comparing the splines to Gaussian processes in Bambi and PyMC

Our interest in these models stems from their flexibility to approximate functions of arbitrary complexity. We’ll see how the methods work in the case of relatively straightforward toy example and then we’ll apply each of the methods to deriving insights into the functional form of insurance loss curves. In this application we adapt a data set discussed in Mick Cooney’s Stan case study to demonstrate the power of hierarchical spline models. Throughout we’ll draw on the discussion of these methods in Osvaldo Martin’s “Bayesian Analysis with Python” for practical details implementing these models.

All of these methods need to be assessed with respect to their in-sample model fit and their out of sample performance. How can we best calibrate the model fits to perform reasonably well out of sample?

Generalised Additive models

The canonical reference for GAMs is Simon Wood’s “Generalised Additive Models: An Introduction with R” which outlines in some detail the theoretical background of splines and univariate smoothers. The book stresses the trade-offs between the flexibility of splines and the need for cross-validation and penalised estimation methods for spline based modelling.

In R these penalised models fits can be achieved in mgcv which incorporates a wilkinson like formula syntax for model specification: y ~ s(x) + s(x1). The closest implementation in python is available in PyGam and we will adopt this package to illustrate the application of smoothing terms and the penalities.

PyGAM and Penalised Fits

In the next code block we load in an example data set on which to demonstrate univariate smoothing patterns using penalised splines. These models can be optimised by fitting differing strength penalities over a varying number of splines

X, y = mcycle(return_X_y=True)

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(X, y)
ax.set_ylabel("Acceleration")
ax.set_xlabel("Time Step")
ax.set_title("Crash Test Dummy Acceleration \n Simulated Motorcycle Crash", fontsize=20)

Text(0.5, 1.0, 'Crash Test Dummy Acceleration \n Simulated Motorcycle Crash')

Next we fit a number of different models to account for the herky-jerky of the data generating processs. We vary the parameterisations to see how the numbers of splines and strength of the penalty help account for the variation in \(y\) over the support of \(X\).

gam1 = LinearGAM(s(0, n_splines=5)).fit(X, y)
gam2 = LinearGAM(s(0, n_splines=7)).fit(X, y)
gam3 = LinearGAM(s(0, n_splines=10)).fit(X, y)
gam4 = LinearGAM(s(0, n_splines=15)).fit(X, y)
gam5 = LinearGAM(s(0, lam=.1)).fit(X, y)
gam6 = LinearGAM(s(0, lam=.5)).fit(X, y)
gam7 = LinearGAM(s(0, lam=5)).fit(X, y)
gam8 = LinearGAM(s(0, lam=15)).fit(X, y)


def plot_fit(gam, X, y, ax, t, c1='b', c2='r'):
    XX = gam.generate_X_grid(term=0, n=500)

    ax.plot(XX, gam.predict(XX), color=c2, linestyle='--')
    ax.plot(XX, gam.prediction_intervals(XX, width=.95), color=c1, ls='--')

    ax.scatter(X, y, facecolor='gray', edgecolors='none')
    ax.set_title(f"""95% prediction interval with {t} \n LL: {gam.statistics_['loglikelihood']}""");

fig, axs = plt.subplots(4,2, figsize=(10, 20))
axs = axs.flatten()
titles = ['5_splines', '7_splines', '10_splines', '15_splines',
'lam=.1', 'lam=.5', 'lam=5', 'lam=15']
gs = [gam1, gam2, gam3, gam4, gam5, gam6, gam7, gam8]
for ax, g, t in zip(axs, gs, titles):
    plot_fit(g, X, y, ax, t)

Here we’ve seen the PyGAM package applied to fitting our model to the data. In the formula specification we see y ~ s(i) where i denotes the index of the column variable in the X data.

Over the range of of the x-axis we can see how the conditional expectation is more or less well fit to the data depending on how the penalities and complexity of the model is specified.

Optimising The Parameter Setting

We can see from the model summary what is going on under the hood. For a given model specification the summary will report a number of model-fit statistics such as the log-likelihood and the AIC.

## Naive Model manually specified splines
gam_raw = LinearGAM(s(0)).fit(X, y)
print("log_likelihood:", gam_raw.statistics_['loglikelihood'])
print("AIC:", gam_raw.statistics_['AIC'])

log_likelihood: -954.5558010528869
AIC: 1931.5965064139743

The question then becomes what changes are induced in the model as we seek to optimise these model fit statistics.

## model optimised
gam = LinearGAM(s(0),  fit_intercept=False)
gam.gridsearch(X, y)
gam.summary()

  0% (0 of 11) |                         | Elapsed Time: 0:00:00 ETA:  --:--:--

  9% (1 of 11) |##                       | Elapsed Time: 0:00:00 ETA:   0:00:00

100% (11 of 11) |########################| Elapsed Time: 0:00:00 Time:  0:00:00

LinearGAM                                                                                                 
=============================================== ==========================================================
Distribution:                        NormalDist Effective DoF:                                     11.8135
Link Function:                     IdentityLink Log Likelihood:                                  -952.2409
Number of Samples:                          133 AIC:                                             1930.1088
                                                AICc:                                            1933.0789
                                                GCV:                                              609.3811
                                                Scale:                                            512.7965
                                                Pseudo R-Squared:                                   0.7984
==========================================================================================================
Feature Function                  Lambda               Rank         EDoF         P > x        Sig. Code   
================================= ==================== ============ ============ ============ ============
s(0)                              [0.2512]             20           11.8         1.11e-16     ***         
==========================================================================================================
Significance codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

WARNING: Fitting splines and a linear function to a feature introduces a model identifiability problem
         which can cause p-values to appear significant when they are not.

WARNING: p-values calculated in this manner behave correctly for un-penalized models or models with
         known smoothing parameters, but when smoothing parameters have been estimated, the p-values
         are typically lower than they should be, meaning that the tests reject the null too readily.

Fortunately, this routine can be performed directly and results in the following differences between the naive and optimised model.

Plot the GAM fits

fig, ax = plt.subplots(figsize=(10, 6))
plot_fit(gam_raw, X, y, ax, "Unoptimised Fit", c1='orange', c2='green')
plot_fit(gam, X, y, ax, "Optimised Fit")

This is all well and good. We’ve seen an approach to modelling that can capture eccentric patterns in raw data. But how does it work and why should we care?

If you’re familiar with the nomenclature of machine learning, you should think of spline modelling as a variety of feature creation. It generates “synthetic” features over the range of the observed variable. These synthetic features are the splines in question.

Digression on the usage of “Splines”

The history of the term “spline” is related to the history of draftmanship. Historically splines were thin strips of flexible wood or plastic that could be bent or shaped around a weight or “knot” points to express a traceable curve over the space of a “numberline”. The elastic nature of the spline material allowed it to be bent around the knot points of curvature expressing a smooth or continuous bend.

The mathematical technique apes these properties by defining a curve over in an analogous way. We specify “knots” to carve up the support of the random variable \(X\) into portions that require different weighting schemes to represent the outcome \(y\)

Extracting the Splines

We can extract the spline features used in the PyGAM by invoking the following commands to first identify the knot points and create the b-spline basis appropriate for the variable \(X\).

knot_edges=utils.gen_edge_knots(X,dtype='numerical')
knots=np.linspace(knot_edges[0],knot_edges[-1],len(gam.coef_))

splines = utils.b_spline_basis(X, edge_knots=knot_edges, sparse=False)

splines_df = pd.DataFrame(splines, columns=[f'basis_{i}' for i in range(len(gam.coef_))])

splines_df.head(10)

	basis_0	basis_1	basis_2	basis_3	basis_4	basis_5	basis_6	basis_7	basis_8	basis_9	basis_10	basis_11	basis_12	basis_13	basis_14	basis_15	basis_16	basis_17	basis_18	basis_19
0	0.166667	0.666667	0.166667	0.000000	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
1	0.132997	0.661606	0.205334	0.000063	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
2	0.059688	0.594827	0.341426	0.004059	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
3	0.030095	0.518726	0.437481	0.013698	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
4	0.012374	0.428013	0.527144	0.032470	0.000000	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
5	0.000000	0.040337	0.551431	0.399315	0.008917	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
6	0.000000	0.018232	0.465467	0.492630	0.023671	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
7	0.000000	0.011137	0.418489	0.535407	0.034967	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
8	0.000000	0.000014	0.189310	0.664817	0.145859	0.000000	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
9	0.000000	0.000000	0.120914	0.656897	0.222015	0.000174	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

These spline features range the extent of the of covariate space \(X\) defining “partitions” of the space. The model “learns” to capture the shape of the outcome variable \(y\) by figuring out how to weight the different portion of this spline basis matrix i.e. the linear combination of this basis matrix with the derived coefficients is a model of our outcome variable.

Note how each row is 0 everywhere except within the columns that represent a partition of \(X\). Additionally each row sum to unity. These properties are important because they ensure that any weighted combination of this basis represents the outcome variable in a controlled and quite granular manner. The more splines we use the more control we have of the representation.

Plotting the Weighted Spline

ax = splines_df.dot(gam.coef_).plot(title='Weighted Splines', label='Weighted Combination of Spline Basis', figsize=(10, 6))
ax.set_ylabel("Acceleration")
ax.set_xlabel("Time Steps")
ax.legend();

In this manner we can see how the specification of a spline basis can help us model eccentric curves in an outcome space.

Next we’ll see how to more directly work with the specification of basis splines, passing these feature matrices into Bambi models.

Bayesian Splines with bambi

Under the hood Bambi makes use of the patsy package formula syntax to specify spline basis terms.

knots_6 = np.linspace(0, np.max(X), 6+2)[1:-1]
knots_10 = np.linspace(0, np.max(X), 10+2)[1:-1]
knots_15 = np.linspace(0, np.max(X), 15+2)[1:-1]

df = pd.DataFrame({'X': X.flatten(), 'y': y})
formula1 = 'bs(X, degree=0, knots=knots_6)'
formula2 = 'bs(X, degree=1, knots=knots_6, intercept=True)'
formula3 = 'bs(X, degree=3, knots=knots_6, intercept=True)'
formula4 = 'bs(X, degree=3, knots=knots_10, intercept=True)'
formula5 = 'bs(X, degree=3, knots=knots_15, intercept=True)'
model_spline1 = bmb.Model(f"y ~ {formula1}", df)
model_spline2 = bmb.Model(f"y ~ {formula2}", df)
model_spline3 = bmb.Model(f"y ~ {formula3}", df)
model_spline4 = bmb.Model(f"y ~ {formula4}", df)
model_spline5 = bmb.Model(f"y ~ {formula5}", df)
model_spline5

       Formula: y ~ bs(X, degree=3, knots=knots_15, intercept=True)
        Family: gaussian
          Link: mu = identity
  Observations: 133
        Priors: 
    target = mu
        Common-level effects
            Intercept ~ Normal(mu: -25.5459, sigma: 226.5005)
            bs(X, degree=3, knots=knots_15, intercept=True) ~ Normal(mu: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
                0. 0. 0. 0. 0. 0. 0. 0.], sigma: [1209.5913 1309.9059 1516.6254  914.7421  803.2567
                526.1027  542.4563
              739.4915  590.0538  692.802   786.9096  707.6947  852.6529  792.6311
             1046.3175 1246.4855 1367.9865 1595.1739 1391.2518])
        
        Auxiliary parameters
            sigma ~ HalfStudentT(nu: 4.0, sigma: 48.14)

As you can see here we are specifying a range of different degrees of spline basis. The different degrees corresspond to the smoothness of the overlapping splines. The degree=0 splines mean we specify a piecewise constant basis i.e. 0 or 1 within each region of the partition. But we can add more degrees to see more flexible representations of the space.

model_spline5.build()
model_spline5.graph()

This will become clearer as we plot the various spline basis matrices below.

Plot the Spline Basis

The below functions extract the basis specification from each model and plots the basis design for an increasingly complex series of spline basis matrices.

def plot_spline_basis(basis, X, ax, title="Spline Basis"):
    df = (
        pd.DataFrame(basis)
        .assign(X=X)
        .melt("X", var_name="basis_idx", value_name="y")
    )


    for idx in df.basis_idx.unique():
        d = df[df.basis_idx == idx]
        ax.plot(d["X"], d["y"])
    
    ax.set_title(title)
    return ax

def plot_knots(knots, ax):
    for knot in knots:
        ax.axvline(knot, color="0.1", alpha=0.4)
    return ax



fig, axs = plt.subplots(5, 1, figsize=(9, 20))
axs = axs.flatten()
axs.flatten()
B1 = model_spline1.response_component.design.common[formula1]
plot_spline_basis(B1, df["X"].values, ax=axs[0], title="Piecewise Constant Basis")
plot_knots(knots_6, axs[0]);

B2 = model_spline2.response_component.design.common[formula2]
ax = plot_spline_basis(B2, df["X"].values, axs[1], 
title="Piecewise Linear Basis")
plot_knots(knots_6, axs[1]);

B3 = model_spline3.response_component.design.common[formula3]
ax = plot_spline_basis(B3, df["X"].values, axs[2], 
title="Cubic Spline Basis (6 Knots)")
plot_knots(knots_6, axs[2]);

B4 = model_spline4.response_component.design.common[formula4]
ax = plot_spline_basis(B4, df["X"].values, axs[3], 
title="Cubic Spline Basis (10 Knots)")
plot_knots(knots_10, axs[3]);


B5 = model_spline5.response_component.design.common[formula5]
ax = plot_spline_basis(B5, df["X"].values, axs[4], 
title="Cubic Spline Basis (15 Knots)")
plot_knots(knots_15, axs[4]);

Now that we’ve seen the nature of the modelling splines we’ll fit each model to the data and plot how the weighted splines matrices are able to represent the raw data.

Fit the Individual Spline Models

idata_spline1 = model_spline1.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

idata_spline2 = model_spline2.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

idata_spline3 = model_spline3.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

idata_spline4 = model_spline4.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

idata_spline5 = model_spline5.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})

Plotting the fits

Plot the Weighted Mean

def plot_weighted_splines(B, idata, formula, ax, knots):
    posterior_stacked = az.extract(idata)
    wp = posterior_stacked[formula].mean("sample").values

    plot_spline_basis(B * wp.T, df["X"].values, ax)
    ax.plot(df.X.values, np.dot(B, wp.T), color="black", lw=3, label='Weighted Splines')
    plot_knots(knots, ax);
    ax.legend()



fig, axs = plt.subplots(5, 1, figsize=(10, 20))
axs = axs.flatten()
axs.flatten()

plot_weighted_splines(B1, idata_spline1, formula1, axs[0], knots_6)
plot_weighted_splines(B2, idata_spline2, formula2, axs[1], knots_6)
plot_weighted_splines(B3, idata_spline3, formula3, axs[2], knots_6)
plot_weighted_splines(B4, idata_spline4, formula4, axs[3], knots_10)
plot_weighted_splines(B5, idata_spline5, formula5, axs[4], knots_15)

Here we can see how the models with increasingly complex splines are more exactly able to fit the herky jerky trajectory of the outcome variable in each interval. The fewer the intervals, the less flexibility available to the model.

Compare Model Fits

As before we can evaluate these model fits and compare them based on leave-one-out cross validation scores and information theoretic complexity measures.

fig, axs = plt.subplots(2, 1, figsize=(9, 16))
axs = axs.flatten()
models_dict = {"piecewise_constant": idata_spline1, "piecewise_linear": idata_spline2, "cubic_bspline": idata_spline3, "cubic_bspline_10": idata_spline4, 
"cubic_bspline_15": idata_spline5}
df_compare = az.compare(models_dict)
az.plot_compare(df_compare, ax=axs[0])
az.plot_compare(az.compare(models_dict, 'waic'), ax=axs[1])

df_compare

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
cubic_bspline_10	0	-612.274608	11.462200	0.000000	9.030402e-01	9.689854	0.000000	True	log
cubic_bspline_15	1	-616.252712	15.566839	3.978104	2.408230e-14	9.833120	1.374420	True	log
cubic_bspline	2	-634.729956	8.698337	22.455348	5.894750e-03	8.897139	6.637243	True	log
piecewise_constant	3	-643.781042	6.981098	31.506433	5.164397e-02	9.770740	8.719803	False	log
piecewise_linear	4	-647.267410	6.170543	34.992802	3.942108e-02	7.902853	8.094122	False	log

Here we see that the extra complexity of using 15 splines leads to slightly worse performance measures than the less complex but seemingly adequate 10 splines.

new_data = pd.DataFrame({"X": np.linspace(df.X.min(), df.X.max(), num=500)})
    
model_spline4.predict(idata_spline4, data=new_data, 
kind='pps', inplace=True)

idata_spline4

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:                                              (chain: 4, draw: 1000,
                                                            bs(X, degree=3, knots=knots_10, intercept=True)_dim: 14,
                                                            y_obs: 500)
  Coordinates:
    * chain                                                (chain) int64 0 1 2 3
    * draw                                                 (draw) int64 0 ... 999
    * bs(X, degree=3, knots=knots_10, intercept=True)_dim  (bs(X, degree=3, knots=knots_10, intercept=True)_dim) int64 ...
    * y_obs                                                (y_obs) int64 0 ... 499
  Data variables:
      Intercept                                            (chain, draw) float64 ...
      bs(X, degree=3, knots=knots_10, intercept=True)      (chain, draw, bs(X, degree=3, knots=knots_10, intercept=True)_dim) float64 ...
      y_sigma                                              (chain, draw) float64 ...
      y_mean                                               (chain, draw, y_obs) float64 ...
  Attributes:
      created_at:                  2024-04-07T11:51:54.606170
      arviz_version:               0.17.0
      inference_library:           pymc
      inference_library_version:   5.10.3
      sampling_time:               7.8021931648254395
      tuning_steps:                1000
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • bs(X, degree=3, knots=knots_10, intercept=True)_dim: 14
    • y_obs: 500
  • Coordinates: (4)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • bs(X, degree=3, knots=knots_10, intercept=True)_dim
      (bs(X, degree=3, knots=knots_10, intercept=True)_dim)
      int64
      0 1 2 3 4 5 6 7 8 9 10 11 12 13

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

  • Data variables: (4)
    • Intercept
      (chain, draw)
      float64
      120.7 -78.41 ... -232.5 -643.4

      array([[ 120.72584539,  -78.40580422, -118.60611945, ...,  121.36419414,
                52.12206675,   90.06865862],
             [ 291.01103275,  312.46955149,  308.2759726 , ..., -220.74994434,
              -221.19586572,  -78.03959137],
             [ -57.96616669,  -65.14951557,  -64.54323192, ..., -482.63520675,
              -482.65677386, -573.52625484],
             [  20.74229491,   77.83153878, -247.93379991, ...,   65.82607356,
              -232.50087738, -643.35960726]])

    • bs(X, degree=3, knots=knots_10, intercept=True)
      (chain, draw, bs(X, degree=3, knots=knots_10, intercept=True)_dim)
      float64
      -122.5 -104.7 ... 637.6 658.6

      array([[[-122.51456067, -104.70757172, -175.78361558, ...,
               -174.58434243,  -97.34264578, -102.47263963],
              [ 105.50470928,  124.7220653 ,   33.02663892, ...,
                 89.63888488,   71.98876682,   99.43913998],
              [ 173.8138609 ,  122.93634145,   60.63548382, ...,
                126.34767578,  125.72172157,  124.27519687],
              ...,
              [-129.48947665, -133.8751114 ,  -88.17353256, ...,
               -122.90446001,  -98.68979663, -149.89281942],
              [ -62.46559992,  -51.11802837,  -95.27067198, ...,
                -41.59539533,  -95.55061917,  -26.39452763],
              [-105.76964538,  -67.0628142 , -134.27233257, ...,
                -49.46032459, -109.82897323,  -90.13466357]],
      
             [[-272.00907129, -272.64520274, -334.1766197 , ...,
               -289.13898071, -296.13606486, -241.3559082 ],
              [-298.72105583, -337.79947478, -316.33303683, ...,
               -321.21665736, -329.86550648, -253.04144959],
              [-323.23183958, -292.06228196, -314.90163154, ...,
               -304.68777643, -297.05477119, -318.64761093],
      ...
              [ 477.17587705,  499.6239322 ,  461.37436802, ...,
                477.05862047,  478.47905055,  514.22366147],
              [ 452.27949475,  523.25026011,  490.73871756, ...,
                448.47506477,  481.63415724,  495.25484626],
              [ 574.13363756,  543.05517393,  552.35598572, ...,
                562.15450166,  571.92641053,  574.61599379]],
      
             [[   7.15842977,  -35.50085456,  -39.01712756, ...,
                  2.39622912,  -52.60133213,   -2.64882026],
              [ -93.19055724,  -69.60274766,  -78.62672525, ...,
                -63.09693048,  -77.05374075,  -50.27036963],
              [ 239.28913699,  274.1687879 ,  231.51398046, ...,
                261.3116388 ,  273.46650709,  260.73531636],
              ...,
              [ -34.3785626 ,  -80.18013599,  -55.9986712 , ...,
                -49.31585595,  -47.86184608,  -51.1751459 ],
              [ 211.8392298 ,  231.02249938,  158.52862855, ...,
                230.46474541,  226.78091415,  233.24213068],
              [ 641.23615142,  625.0061391 ,  667.13099333, ...,
                642.96378814,  637.6410361 ,  658.56934271]]])

    • y_sigma
      (chain, draw)
      float64
      24.99 21.28 22.98 ... 25.23 21.84

      array([[24.98537118, 21.28410029, 22.98278315, ..., 22.15024605,
              25.83916891, 24.71910159],
             [22.65742272, 22.71266799, 22.69326195, ..., 20.55187844,
              20.58224527, 22.16901817],
             [21.74245967, 21.66477858, 21.8756041 , ..., 22.68001903,
              22.65695398, 21.51621838],
             [21.61516254, 22.25846889, 21.11554604, ..., 21.2788373 ,
              25.23339084, 21.84058935]])

    • y_mean
      (chain, draw, y_obs)
      float64
      -1.789 0.1031 1.626 ... 13.91 15.21

      array([[[-1.78871528e+00,  1.03085485e-01,  1.62569616e+00, ...,
                1.86747780e+01,  1.85204557e+01,  1.82532058e+01],
              [ 2.70989051e+01,  2.91166817e+01,  3.06889981e+01, ...,
                1.77447137e+01,  1.93417554e+01,  2.10333358e+01],
              [ 5.52077414e+01,  4.93863093e+01,  4.38241724e+01, ...,
                5.84613354e+00,  5.75920763e+00,  5.66907742e+00],
              ...,
              [-8.12528250e+00, -8.54685180e+00, -8.79218142e+00, ...,
               -2.23722757e+01, -2.53676528e+01, -2.85286253e+01],
              [-1.03435332e+01, -9.13608318e+00, -8.15986093e+00, ...,
                1.74678855e+01,  2.14722713e+01,  2.57275391e+01],
              [-1.57009868e+01, -1.14520878e+01, -7.74519669e+00, ...,
               -2.30115247e+00, -1.24812891e+00, -6.60049429e-02]],
      
             [[ 1.90019615e+01,  1.88337273e+01,  1.84856681e+01, ...,
                4.30189011e+01,  4.62607691e+01,  4.96551245e+01],
              [ 1.37484957e+01,  9.38608034e+00,  5.43201445e+00, ...,
                5.01187900e+01,  5.46671171e+01,  5.94281019e+01],
              [-1.49558670e+01, -1.14852377e+01, -8.35705538e+00, ...,
               -7.76904017e+00, -9.03696042e+00, -1.03716383e+01],
      ...
              [-5.45932970e+00, -2.99400930e+00, -8.40689790e-01, ...,
                2.72430528e+01,  2.93696997e+01,  3.15884547e+01],
              [-3.03772791e+01, -2.24442507e+01, -1.52319734e+01, ...,
                1.08571644e+01,  1.17310518e+01,  1.25980724e+01],
              [ 6.07382720e-01, -2.87409118e+00, -6.05364704e+00, ...,
                7.37590235e-01,  9.16407247e-01,  1.08973896e+00]],
      
             [[ 2.79007247e+01,  2.30962227e+01,  1.86562455e+01, ...,
                1.21692488e+01,  1.50304240e+01,  1.80934746e+01],
              [-1.53590185e+01, -1.27196519e+01, -1.03147655e+01, ...,
                2.43435821e+01,  2.59084898e+01,  2.75611691e+01],
              [-8.64466291e+00, -4.78781691e+00, -1.36587402e+00, ...,
                1.43155480e+01,  1.35833254e+01,  1.28015165e+01],
              ...,
              [ 3.14475110e+01,  2.63328858e+01,  2.16930058e+01, ...,
                1.50491956e+01,  1.48554961e+01,  1.46509277e+01],
              [-2.06616476e+01, -1.86172004e+01, -1.69580267e+01, ...,
               -3.42008269e-02,  3.42761755e-01,  7.41253307e-01],
              [-2.12345584e+00, -3.88306683e+00, -5.37296237e+00, ...,
                1.26814258e+01,  1.39147124e+01,  1.52097355e+01]]])

  • Indexes: (4)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • bs(X, degree=3, knots=knots_10, intercept=True)_dim
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13], dtype='int64', name='bs(X, degree=3, knots=knots_10, intercept=True)_dim'))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='y_obs', length=500))

  • Attributes: (8)

    created_at :

    2024-04-07T11:51:54.606170

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    sampling_time :

    7.8021931648254395

    tuning_steps :

    1000

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• posterior_predictive

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, y_obs: 500)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 492 493 494 495 496 497 498 499
  Data variables:
      y        (chain, draw, y_obs) float64 15.59 -4.381 47.53 ... -4.079 35.94
  Attributes:
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_obs: 500
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

  • Data variables: (1)
    • y
      (chain, draw, y_obs)
      float64
      15.59 -4.381 47.53 ... -4.079 35.94

      array([[[ 1.55911053e+01, -4.38131720e+00,  4.75347119e+01, ...,
                6.97130908e+00,  6.29469988e+01, -2.07532892e+00],
              [ 5.49023205e+01,  4.25591202e+01,  2.77278709e+01, ...,
                3.18060001e+01,  3.17401064e+01,  3.98503541e+01],
              [ 6.90998361e+01,  5.30825171e+01,  1.18360171e+01, ...,
                3.80806793e+01,  2.80825698e+01, -1.94595858e+01],
              ...,
              [-1.59396996e+01, -8.22594303e+00,  7.39448109e-01, ...,
               -5.39676247e+01, -2.65960195e+01, -5.73189493e+00],
              [-1.80704451e+01,  4.08792358e-02, -1.62237628e+01, ...,
                8.00845640e-01,  1.05344423e+02, -9.54887552e-01],
              [-4.11437263e+01,  1.33121241e+01, -1.63066215e+01, ...,
               -4.46076712e+00, -2.25143195e+01,  4.07122304e+01]],
      
             [[ 4.55239057e+01,  1.47871864e+01, -1.07407809e+00, ...,
                2.76735477e+01,  3.01199413e+01,  6.46266446e+01],
              [ 2.76821655e+01,  1.31783880e+01, -3.87877477e+00, ...,
                1.53685075e+01,  5.31560080e+01,  5.94708021e+00],
              [-3.08297987e+01,  2.35201095e+00, -3.45838051e+00, ...,
                2.77489039e+01, -1.66144543e+01,  4.05676592e+00],
      ...
              [-1.63989040e+01, -2.20701423e+01,  2.92735689e+01, ...,
                4.24685161e+01,  5.07346598e+01,  7.04563277e+00],
              [ 2.15423076e+01, -6.96107085e+01, -1.48309782e+01, ...,
                7.86635105e+00,  5.82012257e+00,  1.71137207e+01],
              [-4.21397140e+01, -1.57905523e+01, -2.41613399e+01, ...,
               -3.41975384e+01, -3.24626786e+00,  4.46172972e+01]],
      
             [[ 1.23310003e+01,  4.57747692e+01,  1.91685659e+01, ...,
                3.76472751e+01,  2.07141916e+01,  2.10898490e+01],
              [-3.12560553e+01,  6.12037270e+00, -3.60576274e+01, ...,
                5.40941270e+01,  3.84413068e+01, -4.20487482e+00],
              [-1.47493290e+01, -2.02306624e+01,  2.12970017e+01, ...,
                1.71821453e+01,  6.90235209e+00, -2.22985944e+00],
              ...,
              [ 7.85583846e+00,  4.32391032e+00,  4.10602234e+01, ...,
                1.79621949e+01,  2.73837203e+01,  2.79651899e+01],
              [-1.60331950e+01, -1.45628176e+01,  4.29011105e+01, ...,
                4.80992916e+01, -4.24058830e+01, -2.47372278e+01],
              [-2.18507540e+01, -3.73417596e+01,  1.66484835e+01, ...,
                4.89606803e+01, -4.07926338e+00,  3.59437020e+01]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='y_obs', length=500))

  • Attributes: (2)

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• log_likelihood

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, y_obs: 133)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 125 126 127 128 129 130 131 132
  Data variables:
      y        (chain, draw, y_obs) float64 -4.14 -4.143 -4.178 ... -4.008 -4.024
  Attributes:
      created_at:                  2024-04-07T11:51:54.709906
      arviz_version:               0.17.0
      inference_library:           pymc
      inference_library_version:   5.10.3
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_obs: 133
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 128 129 130 131 132

      array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
              14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
              28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
              42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
              56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
              70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
              84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
              98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
             112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
             126, 127, 128, 129, 130, 131, 132])

  • Data variables: (1)
    • y
      (chain, draw, y_obs)
      float64
      -4.14 -4.143 ... -4.008 -4.024

      array([[[-4.13979163, -4.14290023, -4.17769349, ..., -4.18873926,
               -4.2140095 , -4.18292323],
              [-4.7874173 , -5.08745265, -5.35182326, ..., -4.07333741,
               -4.00381829, -4.09475182],
              [-6.93880952, -6.07172295, -4.48896904, ..., -4.06944286,
               -4.13680171, -4.07764243],
              ...,
              [-4.08406773, -4.07347472, -4.04134234, ..., -4.01843447,
               -4.12572447, -5.58505078],
              [-4.25095186, -4.20784082, -4.18040134, ..., -4.89211473,
               -4.37232929, -4.33994759],
              [-4.32823927, -4.16795339, -4.16177912, ..., -4.26197164,
               -4.1286758 , -4.22135968]],
      
             [[-4.39110458, -4.42381059, -4.31733995, ..., -4.09077766,
               -4.12720544, -5.51743881],
              [-4.22506921, -4.09580767, -4.08371844, ..., -4.2814755 ,
               -4.04651359, -6.34326444],
              [-4.25817634, -4.09755781, -4.07284227, ..., -4.0645694 ,
               -4.09961529, -4.47210136],
      ...
              [-4.06939374, -4.04042752, -4.1131353 , ..., -4.14023954,
               -4.07039226, -4.46454973],
              [-4.93820875, -4.26629647, -4.24927114, ..., -4.36538217,
               -4.04350644, -4.04291443],
              [-3.98814397, -4.00648981, -4.24790848, ..., -4.21120846,
               -3.98780992, -4.08749451]],
      
             [[-4.8254083 , -4.45449574, -4.00572678, ..., -4.49505401,
               -4.06764361, -4.0508328 ],
              [-4.25973202, -4.11195429, -4.02401481, ..., -4.03595185,
               -4.14220154, -4.30857712],
              [-4.05275145, -3.96948228, -4.15676226, ..., -4.05655505,
               -4.54214391, -3.97390066],
              ...,
              [-5.06871224, -4.60468571, -4.03137391, ..., -4.01366979,
               -4.38891768, -3.99388898],
              [-4.48234116, -4.34682584, -4.24183716, ..., -4.32301611,
               -4.1482937 , -4.22498708],
              [-4.00743504, -4.01790161, -4.04838083, ..., -4.16163376,
               -4.00754326, -4.02402652]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             123, 124, 125, 126, 127, 128, 129, 130, 131, 132],
            dtype='int64', name='y_obs', length=133))

  • Attributes: (6)

    created_at :

    2024-04-07T11:51:54.709906

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:                (chain: 4, draw: 1000)
  Coordinates:
    * chain                  (chain) int64 0 1 2 3
    * draw                   (draw) int64 0 1 2 3 4 5 ... 994 995 996 997 998 999
  Data variables: (12/17)
      lp                     (chain, draw) float64 -719.5 -721.5 ... -722.1 -722.2
      index_in_trajectory    (chain, draw) int64 -11 -49 18 4 23 ... -6 -10 57 73
      acceptance_rate        (chain, draw) float64 0.9972 0.6356 ... 0.2066 0.8484
      energy                 (chain, draw) float64 728.2 728.6 ... 729.9 729.2
      step_size              (chain, draw) float64 0.03411 0.03411 ... 0.04662
      n_steps                (chain, draw) float64 47.0 119.0 55.0 ... 127.0 191.0
      ...                     ...
      tree_depth             (chain, draw) int64 6 7 6 4 6 6 7 5 ... 5 5 3 4 5 7 8
      perf_counter_start     (chain, draw) float64 1.549e+06 ... 1.549e+06
      perf_counter_diff      (chain, draw) float64 0.002158 0.005146 ... 0.008517
      energy_error           (chain, draw) float64 -1.295 0.4866 ... 2.466 -2.027
      diverging              (chain, draw) bool False False False ... False False
      smallest_eigval        (chain, draw) float64 nan nan nan nan ... nan nan nan
  Attributes:
      created_at:                  2024-04-07T11:51:54.612916
      arviz_version:               0.17.0
      inference_library:           pymc
      inference_library_version:   5.10.3
      sampling_time:               7.8021931648254395
      tuning_steps:                1000
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

  • Data variables: (17)
    • lp
      (chain, draw)
      float64
      -719.5 -721.5 ... -722.1 -722.2

      array([[-719.50369554, -721.47690958, -722.19648739, ..., -720.90699258,
              -722.11737354, -719.94485737],
             [-721.38353904, -725.29766925, -716.15058152, ..., -720.66343524,
              -721.87096626, -720.99884138],
             [-716.80750964, -720.48567373, -724.60478714, ..., -722.4035202 ,
              -725.60178312, -723.70808827],
             [-716.6405226 , -719.34578606, -717.52389822, ..., -718.05943729,
              -722.11823193, -722.16079379]])

    • index_in_trajectory
      (chain, draw)
      int64
      -11 -49 18 4 23 ... -4 -6 -10 57 73

      array([[-11, -49,  18, ...,  -1, -41,  24],
             [  1, -12,  -6, ...,  -5,   6,  36],
             [ 10,  -9,   6, ..., -50,  -2, -39],
             [ 23,  20, -79, ..., -10,  57,  73]])

    • acceptance_rate
      (chain, draw)
      float64
      0.9972 0.6356 ... 0.2066 0.8484

      array([[0.9971757 , 0.6355526 , 0.96612276, ..., 0.64747254, 0.71357657,
              0.7219741 ],
             [1.        , 0.27075936, 0.74565784, ..., 0.91641409, 0.47375296,
              0.92149302],
             [0.77945693, 0.83027947, 0.61564097, ..., 0.97194468, 0.83515637,
              0.94039973],
             [0.38984477, 0.64599971, 0.95142926, ..., 0.98129117, 0.20662616,
              0.84836349]])

    • energy
      (chain, draw)
      float64
      728.2 728.6 728.7 ... 729.9 729.2

      array([[728.18596723, 728.63910955, 728.68291867, ..., 726.15946692,
              728.15252462, 728.32755882],
             [727.07623419, 734.70606106, 726.35671505, ..., 729.59093046,
              727.79530258, 731.22971686],
             [725.49576441, 726.18720915, 728.13394683, ..., 733.03357126,
              732.76920572, 733.40052227],
             [724.61509281, 726.47676375, 724.5785299 , ..., 724.49490896,
              729.90111393, 729.17952433]])

    • step_size
      (chain, draw)
      float64
      0.03411 0.03411 ... 0.04662 0.04662

      array([[0.0341106 , 0.0341106 , 0.0341106 , ..., 0.0341106 , 0.0341106 ,
              0.0341106 ],
             [0.0362244 , 0.0362244 , 0.0362244 , ..., 0.0362244 , 0.0362244 ,
              0.0362244 ],
             [0.05352544, 0.05352544, 0.05352544, ..., 0.05352544, 0.05352544,
              0.05352544],
             [0.04661661, 0.04661661, 0.04661661, ..., 0.04661661, 0.04661661,
              0.04661661]])

    • n_steps
      (chain, draw)
      float64
      47.0 119.0 55.0 ... 127.0 191.0

      array([[ 47., 119.,  55., ...,  15.,  63.,  63.],
             [  7.,  31.,  15., ...,  15.,  15., 191.],
             [ 39.,  31.,   7., ...,  63.,   7.,  63.],
             [ 47.,  47., 191., ...,  23., 127., 191.]])

    • process_time_diff
      (chain, draw)
      float64
      0.002157 0.005147 ... 0.008517

      array([[0.002157, 0.005147, 0.002465, ..., 0.000732, 0.002836, 0.002888],
             [0.000388, 0.001422, 0.0007  , ..., 0.000727, 0.000719, 0.008406],
             [0.001798, 0.001443, 0.000367, ..., 0.002697, 0.000363, 0.002839],
             [0.002107, 0.002093, 0.008474, ..., 0.001107, 0.005641, 0.008517]])

    • max_energy_error
      (chain, draw)
      float64
      -1.295 1.177 ... 4.799 -2.488

      array([[-1.29520197,  1.17740215, -0.79436509, ...,  0.83609631,
               1.24649329,  1.51146204],
             [-2.00700873,  4.27155897, -3.63691556, ..., -1.13411164,
               1.77054196,  0.40193416],
             [ 0.86450083,  0.45763186,  1.02444297, ..., -0.56746775,
               0.46006735, -0.90094185],
             [ 3.08329139,  1.02752781,  0.40761122, ..., -0.79104795,
               4.79867847, -2.48802427]])

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • step_size_bar
      (chain, draw)
      float64
      0.04864 0.04864 ... 0.05827 0.05827

      array([[0.04863588, 0.04863588, 0.04863588, ..., 0.04863588, 0.04863588,
              0.04863588],
             [0.04742749, 0.04742749, 0.04742749, ..., 0.04742749, 0.04742749,
              0.04742749],
             [0.04865917, 0.04865917, 0.04865917, ..., 0.04865917, 0.04865917,
              0.04865917],
             [0.05827394, 0.05827394, 0.05827394, ..., 0.05827394, 0.05827394,
              0.05827394]])

    • tree_depth
      (chain, draw)
      int64
      6 7 6 4 6 6 7 5 ... 5 5 5 3 4 5 7 8

      array([[6, 7, 6, ..., 4, 6, 6],
             [3, 5, 4, ..., 4, 4, 8],
             [6, 5, 3, ..., 6, 3, 6],
             [6, 6, 8, ..., 5, 7, 8]])

    • perf_counter_start
      (chain, draw)
      float64
      1.549e+06 1.549e+06 ... 1.549e+06

      array([[1548845.98616321, 1548845.98837463, 1548845.99356979, ...,
              1548849.089309  , 1548849.09008833, 1548849.09297283],
             [1548846.26185408, 1548846.26229283, 1548846.26376208, ...,
              1548849.66763733, 1548849.66841163, 1548849.66917683],
             [1548846.00397346, 1548846.00586817, 1548846.00737038, ...,
              1548849.43576621, 1548849.43851262, 1548849.43892025],
             [1548845.75754225, 1548845.75969971, 1548845.76183979, ...,
              1548848.59581225, 1548848.59696821, 1548848.60265904]])

    • perf_counter_diff
      (chain, draw)
      float64
      0.002158 0.005146 ... 0.008517

      array([[0.0021575 , 0.00514625, 0.00246375, ..., 0.00073146, 0.00283608,
              0.00290662],
             [0.00038792, 0.00142125, 0.00070013, ..., 0.00072817, 0.00071883,
              0.00840567],
             [0.00179892, 0.00144321, 0.00036646, ..., 0.002697  , 0.00036288,
              0.00284875],
             [0.00210708, 0.00209321, 0.00848292, ..., 0.00110692, 0.00564096,
              0.00851658]])

    • energy_error
      (chain, draw)
      float64
      -1.295 0.4866 ... 2.466 -2.027

      array([[-1.29520197,  0.48655253, -0.1923247 , ...,  0.58994059,
               0.44985186, -0.34697224],
             [-1.47586104,  3.20202654, -3.63691556, ..., -1.13411164,
               0.28089823, -0.07197861],
             [-0.40128224,  0.32078424,  0.79654477, ..., -0.16681335,
               0.34947932, -0.33374335],
             [-1.16776277,  0.28870592, -0.14708079, ..., -0.79104795,
               2.46552016, -2.02669323]])

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

  • Attributes: (8)

    created_at :

    2024-04-07T11:51:54.612916

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    sampling_time :

    7.8021931648254395

    tuning_steps :

    1000

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:  (y_obs: 133)
  Coordinates:
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 125 126 127 128 129 130 131 132
  Data variables:
      y        (y_obs) float64 0.0 -1.3 -2.7 0.0 -2.7 ... -2.7 10.7 -2.7 10.7
  Attributes:
      created_at:                  2024-04-07T11:51:54.615175
      arviz_version:               0.17.0
      inference_library:           pymc
      inference_library_version:   5.10.3
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • y_obs: 133
  • Coordinates: (1)
    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 128 129 130 131 132

      array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
              14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
              28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
              42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
              56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
              70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
              84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
              98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
             112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
             126, 127, 128, 129, 130, 131, 132])

  • Data variables: (1)
    • y
      (y_obs)
      float64
      0.0 -1.3 -2.7 ... 10.7 -2.7 10.7

      array([   0. ,   -1.3,   -2.7,    0. ,   -2.7,   -2.7,   -2.7,   -1.3,
               -2.7,   -2.7,   -1.3,   -2.7,   -2.7,   -2.7,   -5.4,   -2.7,
               -5.4,    0. ,   -2.7,   -2.7,    0. ,  -13.3,   -5.4,   -5.4,
               -9.3,  -16. ,  -22.8,   -2.7,  -22.8,  -32.1,  -53.5,  -54.9,
              -40.2,  -21.5,  -21.5,  -50.8,  -42.9,  -26.8,  -21.5,  -50.8,
              -61.7,   -5.4,  -80.4,  -59. ,  -71. ,  -91.1,  -77.7,  -37.5,
              -85.6, -123.1, -101.9,  -99.1, -104.4, -112.5,  -50.8, -123.1,
              -85.6,  -72.3, -127.2, -123.1, -117.9, -134. , -101.9, -108.4,
             -123.1, -123.1, -128.5, -112.5,  -95.1,  -81.8,  -53.5,  -64.4,
              -57.6,  -72.3,  -44.3,  -26.8,   -5.4, -107.1,  -21.5,  -65.6,
              -16. ,  -45.6,  -24.2,    9.5,    4. ,   12. ,  -21.5,   37.5,
               46.9,  -17.4,   36.2,   75. ,    8.1,   54.9,   48.2,   46.9,
               16. ,   45.6,    1.3,   75. ,  -16. ,  -54.9,   69.6,   34.8,
               32.1,  -37.5,   22.8,   46.9,   10.7,    5.4,   -1.3,  -21.5,
              -13.3,   30.8,  -10.7,   29.4,    0. ,  -10.7,   14.7,   -1.3,
                0. ,   10.7,   10.7,  -26.8,  -14.7,  -13.3,    0. ,   10.7,
              -14.7,   -2.7,   10.7,   -2.7,   10.7])

  • Indexes: (1)
    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             123, 124, 125, 126, 127, 128, 129, 130, 131, 132],
            dtype='int64', name='y_obs', length=133))

  • Attributes: (6)

    created_at :

    2024-04-07T11:51:54.615175

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  

Next we plot the posterior predictive distribution of our observed variable and compare against the observed data. Additionally we plot the 89th and 50% HDI.

ax = az.plot_hdi(new_data['X'], idata_spline4['posterior_predictive']['y'], fill_kwargs={'alpha': 0.2, 'color':'firebrick'}, hdi_prob=0.89, figsize=(10, 8))

az.plot_hdi(new_data['X'], idata_spline4['posterior_predictive']['y'], fill_kwargs={'alpha': 0.8, 'color':'firebrick'}, hdi_prob=0.5)

y_mean = idata_spline4['posterior_predictive']['y'].mean(dim=('chain', 'draw'))

ax.plot(new_data['X'], y_mean, label='Expected posterior predictive', color='k')
ax.set_xlabel("Time Point")
ax.set_ylabel("Acceleration")

ax.scatter(df['X'], df['y'], label='Observed Datapoints')

ax.legend()

ax.set_title("Posterior Predictive Distribution \n Based on 10 Knots");

This represents a good a clean model fit to the observed data using univariate spline smoothers. Next we’ll see another alternative approach to model this outcome variable using approximate gaussian processes.

Gaussian processes

The topic of gaussian processes is rich and detailed. Too rich to be fairly covered in this blog post, so we’ll just say that we’re using a method designed for function approximation that makes use of drawing samples from a multivariate normal distribution under a range of different covariance relationships.

These relationships can be somewhat intuitively interrogated by defining different combinations of covariance relationships with priors over the parameters governing the covariance of a sequence of points. For example consider the following parameterisations.

lengthscale = 3
sigma = 13
cov = sigma**2 * pm.gp.cov.ExpQuad(1, lengthscale)

X = np.linspace(0, 60, 200)[:, None]
K = cov(X).eval()

fig, ax = plt.subplots(figsize=(9, 7))
ax.plot(
    X,
    pm.draw(
        pm.MvNormal.dist(mu=np.zeros(len(K)), cov=K, shape=K.shape[0]), draws=10, random_seed=random_seed
    ).T,
)
plt.title(f"Samples from the GP prior \n lengthscale: {3}, sigma: {13}")
plt.ylabel("y")
plt.xlabel("X");

We’ve specified the range of X to reflect the support of the acceleration example and allowed the draws to be informed by a covariance function we have parameterised using the Exponentiated Quadratic kernel:

\[k(x, x') = \mathrm{exp}\left[ -\frac{(x - x')^2}{2 \ell^2} \right]\]

The patterns exhibited show a good range of “wiggliness” that they should be flexible enough to capture the shape of the acceleration, if we can calibrate the posterior of parameters against the observed data.

Priors on Gaussian Priors

Consider the following specification:

fig, ax = plt.subplots(figsize=(9, 7))
ax.hist(pm.draw(pm.InverseGamma.dist(mu=1, sigma=1), 1000), ec='black', bins=30);

fig, ax = plt.subplots(figsize=(9, 7))
ax.hist(pm.draw(pm.Exponential.dist(lam=1), 1000), ec='black', bins=30);

We use these to specify priors on the Hilbert space approximation of gaussian priors available in the Bambi package.

prior_hsgp = {
    "sigma": bmb.Prior("Exponential", lam=1), # amplitude
    "ell": bmb.Prior("InverseGamma", mu=1, sigma=1) # lengthscale
}

# This is the dictionary we pass to Bambi
priors = {
    "hsgp(X, m=10, c=1)": prior_hsgp,
    "sigma": bmb.Prior("HalfNormal", sigma=4)
}
model_hsgp = bmb.Model("y ~ 0 + hsgp(X, m=10, c=1)", df, priors=priors)
model_hsgp

       Formula: y ~ 0 + hsgp(X, m=10, c=1)
        Family: gaussian
          Link: mu = identity
  Observations: 133
        Priors: 
    target = mu
        HSGP contributions
            hsgp(X, m=10, c=1)
                cov: ExpQuad
                sigma ~ Exponential(lam: 1.0)
                ell ~ InverseGamma(mu: 1.0, sigma: 1.0)
        
        Auxiliary parameters
            sigma ~ HalfNormal(sigma: 4.0)

Here we’ve set the m=10 to determine the number of basis vectors used in the Hilbert space approximation. The idea differs in detail from the spline based approximations we’ve seen, but it’s perhaps useful to think of the process in the same vein.

idata = model_hsgp.fit(inference_method="nuts_numpyro",target_accept=0.95, random_seed=121195, 
idata_kwargs={"log_likelihood": True})
print(idata.sample_stats["diverging"].sum().to_numpy())

Compiling...

Compilation time = 0:00:01.082921

Sampling...

Sampling time = 0:00:01.871345

Transforming variables...

Transformation time = 0:00:00.177758

Computing Log Likelihood...

Log Likelihood time = 0:00:00.143160

0

This model fits and the sampling seems to have worked well.

az.plot_trace(idata, backend_kwargs={"layout": "constrained"}, figsize=(9, 15));

The lengthscale and sigma parameters we have learned by calibrating our priors against the data. The degree to which these parameters are meaningful depend a little on how familar you are with covariance matrix kernels and their properties, so we won’t dwell on the point here.

az.summary(idata, var_names=['hsgp(X, m=10, c=1)_ell', 'hsgp(X, m=10, c=1)_sigma', 'y_sigma', 'hsgp(X, m=10, c=1)_weights'])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
hsgp(X, m=10, c=1)_ell	3.301	0.666	2.099	4.567	0.016	0.011	1763.0	1815.0	1.0
hsgp(X, m=10, c=1)_sigma	23.675	2.679	18.925	28.847	0.069	0.049	1497.0	2167.0	1.0
y_sigma	20.589	1.121	18.590	22.751	0.019	0.014	3440.0	2485.0	1.0
hsgp(X, m=10, c=1)_weights[0]	-131.588	13.784	-157.019	-105.189	0.216	0.154	4071.0	3020.0	1.0
hsgp(X, m=10, c=1)_weights[1]	-94.391	18.010	-125.921	-58.020	0.363	0.256	2500.0	2876.0	1.0
hsgp(X, m=10, c=1)_weights[2]	106.121	20.717	69.627	147.247	0.458	0.324	2046.0	2563.0	1.0
hsgp(X, m=10, c=1)_weights[3]	149.289	22.154	106.508	188.107	0.521	0.368	1807.0	2467.0	1.0
hsgp(X, m=10, c=1)_weights[4]	-83.399	22.728	-124.812	-40.406	0.539	0.383	1776.0	2517.0	1.0
hsgp(X, m=10, c=1)_weights[5]	-125.322	23.273	-169.346	-82.187	0.531	0.377	1917.0	2940.0	1.0
hsgp(X, m=10, c=1)_weights[6]	37.835	21.153	-1.712	77.574	0.512	0.362	1713.0	2577.0	1.0
hsgp(X, m=10, c=1)_weights[7]	94.574	20.745	53.167	132.182	0.423	0.299	2393.0	3105.0	1.0
hsgp(X, m=10, c=1)_weights[8]	-1.494	17.184	-31.835	33.063	0.361	0.255	2264.0	2853.0	1.0
hsgp(X, m=10, c=1)_weights[9]	-45.442	15.955	-76.649	-16.636	0.289	0.208	3041.0	3183.0	1.0

But again we can sample from the posterior predictive distribution of the outcome variable

model_hsgp.predict(idata, data=new_data, 
kind='pps', inplace=True)

idata

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:                         (chain: 4, draw: 1000,
                                       hsgp(X, m=10, c=1)_weights_dim: 10,
                                       y_obs: 500)
  Coordinates:
    * chain                           (chain) int64 0 1 2 3
    * draw                            (draw) int64 0 1 2 3 4 ... 996 997 998 999
    * hsgp(X, m=10, c=1)_weights_dim  (hsgp(X, m=10, c=1)_weights_dim) int64 0 ...
    * y_obs                           (y_obs) int64 0 1 2 3 4 ... 496 497 498 499
  Data variables:
      hsgp(X, m=10, c=1)_weights_raw  (chain, draw, hsgp(X, m=10, c=1)_weights_dim) float64 ...
      y_sigma                         (chain, draw) float64 20.86 19.86 ... 21.17
      hsgp(X, m=10, c=1)_sigma        (chain, draw) float64 26.85 24.0 ... 18.06
      hsgp(X, m=10, c=1)_ell          (chain, draw) float64 3.861 3.509 ... 3.76
      hsgp(X, m=10, c=1)_weights      (chain, draw, hsgp(X, m=10, c=1)_weights_dim) float64 ...
      y_mean                          (chain, draw, y_obs) float64 -0.1748 ... ...
      hsgp(X, m=10, c=1)              (chain, draw, y_obs) float64 -0.1748 ... ...
  Attributes:
      created_at:                  2024-04-07T11:52:14.627428
      arviz_version:               0.17.0
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • hsgp(X, m=10, c=1)_weights_dim: 10
    • y_obs: 500
  • Coordinates: (4)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • hsgp(X, m=10, c=1)_weights_dim
      (hsgp(X, m=10, c=1)_weights_dim)
      int64
      0 1 2 3 4 5 6 7 8 9

      array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

  • Data variables: (7)
    • hsgp(X, m=10, c=1)_weights_raw
      (chain, draw, hsgp(X, m=10, c=1)_weights_dim)
      float64
      -1.452 -0.974 ... 0.9659 -0.2453

      array([[[-1.45186697, -0.97401433,  1.51626692, ...,  1.84709969,
                0.92584892, -1.19616826],
              [-1.73799921, -1.55920285,  1.42827017, ...,  1.78787921,
               -1.14159607, -1.66042935],
              [-1.3059213 , -1.22687096,  2.0214052 , ...,  2.80079646,
                0.17462599, -0.06073935],
              ...,
              [-2.26781832, -1.21871881,  1.80249133, ...,  2.65551316,
                0.30040289, -0.7085128 ],
              [-2.23372877, -1.1832454 ,  1.74365597, ...,  2.71376469,
                0.41123845, -0.60671544],
              [-1.96117304, -1.13093178,  2.24600086, ...,  2.3750183 ,
               -0.00901474, -0.75147244]],
      
             [[-1.95999601, -1.27507311,  1.52467611, ...,  3.28813978,
               -0.1598232 , -1.10988783],
              [-1.59314189, -1.44864255,  1.95956808, ...,  2.01431167,
                1.0199955 , -0.8962025 ],
              [-1.8432561 , -1.32415496,  1.80874966, ...,  3.05464118,
                0.81105708, -1.02872434],
      ...
              [-1.94277245, -1.5433581 ,  1.78771407, ...,  2.11464836,
               -0.62470516, -1.01229181],
              [-1.96747626, -1.39055532,  1.24775496, ...,  2.93940091,
                0.46989555, -1.57735663],
              [-2.20654554, -1.042289  ,  1.46402154, ...,  3.37106106,
               -0.11061588, -2.15605928]],
      
             [[-1.68378706, -0.84035677,  1.41475357, ...,  2.46320953,
                0.88793346, -0.56144583],
              [-1.66959934, -1.01014124,  1.6874097 , ...,  2.71030445,
                0.68597924, -0.76198643],
              [-1.54563186, -0.89279988,  1.77115866, ...,  2.3842957 ,
                0.08392828, -0.78058783],
              ...,
              [-2.17148384, -1.59167361,  1.34960308, ...,  2.34802173,
                0.44912308, -1.0698665 ],
              [-1.31376138, -0.63355604,  1.97759132, ...,  1.8841587 ,
               -0.13453184, -1.36605679],
              [-2.0018976 , -1.36470426,  2.15352253, ...,  3.17100177,
                0.96590611, -0.24527631]]])

    • y_sigma
      (chain, draw)
      float64
      20.86 19.86 21.17 ... 20.12 21.17

      array([[20.85613959, 19.86298062, 21.16886643, ..., 19.48940449,
              19.31511941, 19.44304302],
             [18.5908299 , 21.71605063, 21.22255223, ..., 19.44547467,
              19.21042952, 22.0099345 ],
             [20.48425625, 21.17420065, 22.11318822, ..., 18.40825711,
              21.65013089, 22.05150211],
             [21.08143229, 20.7215663 , 22.60954219, ..., 21.15149371,
              20.11622356, 21.16568101]])

    • hsgp(X, m=10, c=1)_sigma
      (chain, draw)
      float64
      26.85 24.0 27.81 ... 26.09 18.06

      array([[26.84572148, 24.00472166, 27.80754218, ..., 23.34108117,
              23.72941745, 22.60241785],
             [22.01917097, 24.17709343, 26.16315265, ..., 22.03189338,
              22.17588023, 24.00832607],
             [24.12239909, 24.36393141, 25.79692715, ..., 24.33866729,
              20.63877848, 22.26540398],
             [25.41092492, 25.05093181, 28.3790263 , ..., 22.30157727,
              26.0901197 , 18.06420147]])

    • hsgp(X, m=10, c=1)_ell
      (chain, draw)
      float64
      3.861 3.509 3.695 ... 3.679 3.76

      array([[3.86057038, 3.50871714, 3.69475981, ..., 3.20961859, 3.07237787,
              2.99194079],
             [4.21093385, 4.11858175, 3.55744155, ..., 3.54737759, 2.43665447,
              3.67009957],
             [3.38790729, 3.56906294, 3.07026931, ..., 3.17401967, 4.77493695,
              4.1757221 ],
             [3.73124202, 3.62463311, 4.06803959, ..., 3.06030257, 3.67928011,
              3.75995625]])

    • hsgp(X, m=10, c=1)_weights
      (chain, draw, hsgp(X, m=10, c=1)_weights_dim)
      float64
      -120.2 -78.54 ... 27.36 -5.933

      array([[[-120.19159663,  -78.54474016,  117.04015633, ...,
                 88.13161831,   38.07210878,  -41.65695052],
              [-122.83658383, -107.83699596,   95.27680851, ...,
                 80.15736607,  -45.26655282,  -57.3944273 ],
              [-109.63228325, -100.55013819,  159.16260856, ...,
                141.94233647,    7.72311493,   -2.30700494],
              ...,
              [-149.23661164,  -78.75788623,  113.01499747, ...,
                119.40081633,   12.18790241,  -25.62639515],
              [-146.28304701,  -76.21158457,  109.23906639, ...,
                125.36332205,   17.28993874,  -22.96014529],
              [-120.75657952,  -68.54675039,  132.60289356, ...,
                105.03375274,   -0.36461133,  -27.50691292]],
      
             [[-138.76216692,  -87.49701007,   99.32062774, ...,
                120.85182044,   -4.92170357,  -28.04716014],
              [-122.53330814, -108.14123925,  139.17924148, ...,
                 82.74939574,   35.37871319,  -25.72818479],
              [-142.94342858, -100.42482592,  132.17632445, ...,
                148.36712045,   34.72137014,  -38.24384896],
      ...
              [-132.58698442, -103.47676242,  116.36868708, ...,
                 99.43921353,  -26.56716678,  -38.47612061],
              [-138.61535955,  -94.11467693,   78.98484564, ...,
                 89.13997713,   11.35087487,  -29.54951483],
              [-157.32953055,  -72.06983662,   96.1818003 , ...,
                126.1515611 ,   -3.47852942,  -55.82170316]],
      
             [[-129.78695588,  -63.20664967,  102.15003445, ...,
                113.4760498 ,   35.60103971,  -19.27405014],
              [-125.10213961,  -73.9587159 ,  118.87369923, ...,
                124.94562328,   27.73901676,  -26.61384578],
              [-138.71497899,  -77.82483202,  147.07295446, ...,
                116.05385128,    3.46342745,  -26.78443203],
              ...,
              [-133.39280566,  -96.17651826,   79.33895456, ...,
                102.02385007,   17.77408526,  -38.1416893 ],
              [-103.26872735,  -48.62806342,  145.87733997, ...,
                 89.78661037,   -5.60101795,  -48.90554932],
              [-110.1016679 ,  -73.21183754,  110.83494684, ...,
                103.40789151,   27.35517734,   -5.93337047]]])

    • y_mean
      (chain, draw, y_obs)
      float64
      -0.1748 -0.6437 ... 4.329e-15

      array([[[-1.74825831e-01, -6.43722945e-01, -1.10198757e+00, ...,
                8.81438724e-01,  4.41420533e-01,  1.00916482e-14],
              [ 3.99710952e+00,  4.36245618e+00,  4.70113794e+00, ...,
               -3.68931837e-01, -1.84461976e-01, -4.21486081e-15],
              [ 1.12371919e+01,  9.85417314e+00,  8.49086985e+00, ...,
               -4.14235314e-01, -2.07469138e-01, -4.74327547e-15],
              ...,
              [ 5.14480643e-01, -3.24968824e-01, -1.15248120e+00, ...,
                2.85647398e-02,  1.43964056e-02,  3.29822057e-16],
              [-1.15486831e+00, -2.03017733e+00, -2.88893220e+00, ...,
               -1.21523155e-01, -6.06829841e-02, -1.38598577e-15],
              [ 1.21777447e+01,  1.16287837e+01,  1.10821348e+01, ...,
                1.95250577e-01,  9.77499077e-02,  2.23450097e-15]],
      
             [[-2.44479507e+00, -2.94292536e+00, -3.43456273e+00, ...,
               -7.42363606e-01, -3.71378424e-01, -8.48736477e-15],
              [-3.36694245e+00, -3.93597519e+00, -4.48819398e+00, ...,
                1.45535040e+00,  7.28337111e-01,  1.66472919e-14],
              [ 1.86823417e-01, -1.03465399e+00, -2.23999952e+00, ...,
                3.24627418e-01,  1.62576930e-01,  3.71683234e-15],
      ...
              [-3.54206467e+00, -3.29770384e+00, -3.05890542e+00, ...,
               -8.07266002e-01, -4.03861987e-01, -9.22984926e-15],
              [-8.28095921e+00, -8.58481681e+00, -8.88022655e+00, ...,
               -1.37004285e-01, -6.83094783e-02, -1.55937910e-15],
              [ 1.05497489e+00,  4.74317876e-01, -1.10703434e-01, ...,
               -4.72148575e-01, -2.35887059e-01, -5.38850915e-15]],
      
             [[-4.02175043e+00, -4.75614076e+00, -5.46523621e+00, ...,
               -1.18443060e-01, -5.89840516e-02, -1.34595499e-15],
              [ 3.62393761e+00,  2.72714782e+00,  1.84762910e+00, ...,
               -9.97541423e-02, -4.96901944e-02, -1.13397909e-15],
              [ 1.18490633e+01,  1.10655715e+01,  1.02860360e+01, ...,
               -5.74333102e-01, -2.87269966e-01, -6.56480568e-15],
              ...,
              [-9.39146019e+00, -9.93167964e+00, -1.04627004e+01, ...,
                1.23581084e-01,  6.21200724e-02,  1.42193352e-15],
              [ 2.48719138e+01,  2.46289484e+01,  2.43707844e+01, ...,
               -9.35532279e-02, -4.65445693e-02, -1.06175818e-15],
              [-5.90472295e-01, -1.42675597e+00, -2.23855095e+00, ...,
                3.78505074e-01,  1.89406640e-01,  4.32904960e-15]]])

    • hsgp(X, m=10, c=1)
      (chain, draw, y_obs)
      float64
      -0.1748 -0.6437 ... 4.329e-15

      array([[[-1.74825831e-01, -6.43722945e-01, -1.10198757e+00, ...,
                8.81438724e-01,  4.41420533e-01,  1.00916482e-14],
              [ 3.99710952e+00,  4.36245618e+00,  4.70113794e+00, ...,
               -3.68931837e-01, -1.84461976e-01, -4.21486081e-15],
              [ 1.12371919e+01,  9.85417314e+00,  8.49086985e+00, ...,
               -4.14235314e-01, -2.07469138e-01, -4.74327547e-15],
              ...,
              [ 5.14480643e-01, -3.24968824e-01, -1.15248120e+00, ...,
                2.85647398e-02,  1.43964056e-02,  3.29822057e-16],
              [-1.15486831e+00, -2.03017733e+00, -2.88893220e+00, ...,
               -1.21523155e-01, -6.06829841e-02, -1.38598577e-15],
              [ 1.21777447e+01,  1.16287837e+01,  1.10821348e+01, ...,
                1.95250577e-01,  9.77499077e-02,  2.23450097e-15]],
      
             [[-2.44479507e+00, -2.94292536e+00, -3.43456273e+00, ...,
               -7.42363606e-01, -3.71378424e-01, -8.48736477e-15],
              [-3.36694245e+00, -3.93597519e+00, -4.48819398e+00, ...,
                1.45535040e+00,  7.28337111e-01,  1.66472919e-14],
              [ 1.86823417e-01, -1.03465399e+00, -2.23999952e+00, ...,
                3.24627418e-01,  1.62576930e-01,  3.71683234e-15],
      ...
              [-3.54206467e+00, -3.29770384e+00, -3.05890542e+00, ...,
               -8.07266002e-01, -4.03861987e-01, -9.22984926e-15],
              [-8.28095921e+00, -8.58481681e+00, -8.88022655e+00, ...,
               -1.37004285e-01, -6.83094783e-02, -1.55937910e-15],
              [ 1.05497489e+00,  4.74317876e-01, -1.10703434e-01, ...,
               -4.72148575e-01, -2.35887059e-01, -5.38850915e-15]],
      
             [[-4.02175043e+00, -4.75614076e+00, -5.46523621e+00, ...,
               -1.18443060e-01, -5.89840516e-02, -1.34595499e-15],
              [ 3.62393761e+00,  2.72714782e+00,  1.84762910e+00, ...,
               -9.97541423e-02, -4.96901944e-02, -1.13397909e-15],
              [ 1.18490633e+01,  1.10655715e+01,  1.02860360e+01, ...,
               -5.74333102e-01, -2.87269966e-01, -6.56480568e-15],
              ...,
              [-9.39146019e+00, -9.93167964e+00, -1.04627004e+01, ...,
                1.23581084e-01,  6.21200724e-02,  1.42193352e-15],
              [ 2.48719138e+01,  2.46289484e+01,  2.43707844e+01, ...,
               -9.35532279e-02, -4.65445693e-02, -1.06175818e-15],
              [-5.90472295e-01, -1.42675597e+00, -2.23855095e+00, ...,
                3.78505074e-01,  1.89406640e-01,  4.32904960e-15]]])

  • Indexes: (4)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • hsgp(X, m=10, c=1)_weights_dim
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='hsgp(X, m=10, c=1)_weights_dim'))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='y_obs', length=500))

  • Attributes: (4)

    created_at :

    2024-04-07T11:52:14.627428

    arviz_version :

    0.17.0

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• posterior_predictive

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, y_obs: 500)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 492 493 494 495 496 497 498 499
  Data variables:
      y        (chain, draw, y_obs) float64 1.602 -43.6 -19.57 ... -0.1009 27.57
  Attributes:
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_obs: 500
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

  • Data variables: (1)
    • y
      (chain, draw, y_obs)
      float64
      1.602 -43.6 ... -0.1009 27.57

      array([[[  1.60168428, -43.59568681, -19.57451424, ...,  61.74364622,
                 0.93530324,  -1.02814658],
              [ 30.04991774,   5.14866514,  15.63647936, ...,   7.79368728,
               -12.71488575, -13.6055619 ],
              [ 11.4762129 ,   8.10095539, -10.73565714, ..., -11.41875954,
                14.39509793,  10.95749635],
              ...,
              [ -3.87077008, -20.74166775,  32.67554916, ...,  10.01886157,
                 3.92077449, -37.4122051 ],
              [ -4.61104543,  12.48885423,  24.52332409, ...,  14.48459552,
                12.68287723,  -6.83099474],
              [ 29.98606836,   0.8475416 ,  23.73368482, ...,   8.25166262,
                 5.7547574 ,  -2.58944686]],
      
             [[ -7.16470775, -46.84884538, -25.46922122, ..., -16.07794441,
                23.94430991,   5.3644025 ],
              [ 31.93450232, -33.22524479, -32.09356185, ...,   3.17930767,
                12.85620708, -11.41908142],
              [ 12.03305536,   8.39450872, -30.30656756, ..., -21.2604651 ,
                39.25314584,  30.69370805],
      ...
              [ 13.60908896, -43.12812545, -25.00944683, ...,  18.16137521,
                 9.37951823,  26.74217056],
              [-13.26058458, -35.16254651,   1.21184678, ...,  -1.76957082,
                31.03973927, -33.84153268],
              [-28.32935144,   6.18624638,  -9.17055657, ...,  22.52832216,
                11.36442558,  25.2035807 ]],
      
             [[ -8.70999277,  14.3345426 ,  30.62329733, ..., -18.70006159,
                -0.1418926 ,  24.07609217],
              [ 26.55579728,  -1.72476756, -21.80701529, ...,  28.8223474 ,
               -26.84141068,  28.91645408],
              [ 26.80641985,  42.08690557,  -4.81274197, ...,   8.18379724,
                 3.68782771, -11.31903714],
              ...,
              [  3.34426769, -11.85502409, -18.08649394, ..., -46.31834546,
               -31.9137845 ,   1.93764836],
              [ 36.48816743,  -1.60509697,  18.9994632 , ...,  10.34894315,
                -5.6655749 ,  -5.5573968 ],
              [ -7.64079423,  10.87234419,   2.15500364, ..., -10.46747553,
                -0.10087789,  27.56732543]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='y_obs', length=500))

  • Attributes: (2)

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• log_likelihood

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, y_obs: 133)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 125 126 127 128 129 130 131 132
  Data variables:
      y        (chain, draw, y_obs) float64 -3.957 -3.957 -3.957 ... -4.013 -4.099
  Attributes:
      created_at:                  2024-04-07T11:52:14.630993
      arviz_version:               0.17.0
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_obs: 133
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 128 129 130 131 132

      array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
              14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
              28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
              42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
              56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
              70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
              84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
              98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
             112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
             126, 127, 128, 129, 130, 131, 132])

  • Data variables: (1)
    • y
      (chain, draw, y_obs)
      float64
      -3.957 -3.957 ... -4.013 -4.099

      array([[[-3.95662204, -3.95668036, -3.95697947, ..., -3.96786393,
               -4.06567925, -4.08819109],
              [-3.92804385, -3.95248402, -4.00438808, ..., -4.19279911,
               -3.90894979, -4.05289002],
              [-4.11236339, -4.08418887, -3.99366348, ..., -4.1948162 ,
               -3.97180293, -4.09921461],
              ...,
              [-3.88915792, -3.88893229, -3.89724674, ..., -4.04248536,
               -3.89847518, -4.03951892],
              [-3.88161419, -3.88254965, -3.9052254 , ..., -4.08668059,
               -3.88211318, -4.03326819],
              [-4.08257196, -4.09265666, -4.04648571, ..., -3.99108035,
               -3.91131912, -4.03785686]],
      
             [[-3.85025382, -3.84763163, -3.85608195, ..., -4.34165517,
               -3.86684608, -4.00723736],
              [-4.0090095 , -4.00707108, -4.017144  , ..., -4.01175855,
               -4.25268066, -4.11837821],
              [-3.97404168, -3.97456237, -4.00764182, ..., -4.0433394 ,
               -4.00511976, -4.10110199],
      ...
              [-3.85025004, -3.83654156, -3.832685  , ..., -4.37637177,
               -3.86484545, -4.00066997],
              [-4.06709934, -4.05434047, -4.05493257, ..., -4.17570336,
               -3.99476606, -4.11607839],
              [-4.01346366, -4.01406212, -4.01257167, ..., -4.30726434,
               -4.01839009, -4.13004222]],
      
             [[-3.98552821, -3.98561292, -4.0079376 , ..., -4.15815157,
               -3.96833515, -4.09613757],
              [-3.96540632, -3.96291088, -3.95021137, ..., -4.13476712,
               -3.95195981, -4.08343265],
              [-4.17463712, -4.17201421, -4.11695573, ..., -4.3197009 ,
               -4.04463544, -4.14929407],
              ...,
              [-4.0692212 , -4.06241625, -4.09034963, ..., -4.09469991,
               -3.98159516, -4.09860353],
              [-4.68482086, -4.73793167, -4.72387693, ..., -4.11913623,
               -3.92229916, -4.06192875],
              [-3.97170872, -3.97200671, -3.98370832, ..., -4.02417539,
               -4.01254991, -4.09910258]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             123, 124, 125, 126, 127, 128, 129, 130, 131, 132],
            dtype='int64', name='y_obs', length=133))

  • Attributes: (4)

    created_at :

    2024-04-07T11:52:14.630993

    arviz_version :

    0.17.0

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:          (chain: 4, draw: 1000)
  Coordinates:
    * chain            (chain) int64 0 1 2 3
    * draw             (draw) int64 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
  Data variables:
      acceptance_rate  (chain, draw) float64 0.9484 0.989 0.9772 ... 0.9905 0.9283
      step_size        (chain, draw) float64 0.1631 0.1631 ... 0.1827 0.1827
      diverging        (chain, draw) bool False False False ... False False False
      energy           (chain, draw) float64 671.4 665.4 671.0 ... 670.8 674.0
      n_steps          (chain, draw) int64 31 31 31 15 15 31 ... 15 31 31 15 15 31
      tree_depth       (chain, draw) int64 5 5 5 4 4 5 5 4 4 ... 4 5 5 4 5 5 4 4 5
      lp               (chain, draw) float64 661.7 661.3 665.7 ... 665.6 668.8
  Attributes:
      created_at:                  2024-04-07T11:52:14.629957
      arviz_version:               0.17.0
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

  • Data variables: (7)
    • acceptance_rate
      (chain, draw)
      float64
      0.9484 0.989 ... 0.9905 0.9283

      array([[0.94837217, 0.98900712, 0.97718268, ..., 0.99599719, 0.97974622,
              0.93727528],
             [0.9824305 , 0.87426419, 0.99008681, ..., 0.99942845, 0.97008375,
              0.99555447],
             [0.99910745, 0.98130743, 0.95550076, ..., 0.97069333, 0.9342073 ,
              0.95033605],
             [0.97363799, 0.97895034, 0.93808177, ..., 0.9514047 , 0.99054107,
              0.92830969]])

    • step_size
      (chain, draw)
      float64
      0.1631 0.1631 ... 0.1827 0.1827

      array([[0.16306409, 0.16306409, 0.16306409, ..., 0.16306409, 0.16306409,
              0.16306409],
             [0.15671706, 0.15671706, 0.15671706, ..., 0.15671706, 0.15671706,
              0.15671706],
             [0.15280013, 0.15280013, 0.15280013, ..., 0.15280013, 0.15280013,
              0.15280013],
             [0.1826574 , 0.1826574 , 0.1826574 , ..., 0.1826574 , 0.1826574 ,
              0.1826574 ]])

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • energy
      (chain, draw)
      float64
      671.4 665.4 671.0 ... 670.8 674.0

      array([[671.39512608, 665.35858794, 670.96295141, ..., 670.72923684,
              668.42111758, 669.24601286],
             [665.04518693, 671.40450254, 672.47055359, ..., 665.8504822 ,
              674.75907256, 672.99992931],
             [672.00375668, 668.16520377, 666.3166769 , ..., 671.79404347,
              671.22701287, 671.74308362],
             [669.46327702, 671.11966704, 667.28825365, ..., 665.11423361,
              670.79366043, 673.9890605 ]])

    • n_steps
      (chain, draw)
      int64
      31 31 31 15 15 ... 31 31 15 15 31

      array([[31, 31, 31, ..., 31, 15, 15],
             [15, 31, 31, ..., 31, 31, 31],
             [31, 31, 15, ..., 15, 31, 15],
             [31, 15, 15, ..., 15, 15, 31]])

    • tree_depth
      (chain, draw)
      int64
      5 5 5 4 4 5 5 4 ... 5 5 4 5 5 4 4 5

      array([[5, 5, 5, ..., 5, 4, 4],
             [4, 5, 5, ..., 5, 5, 5],
             [5, 5, 4, ..., 4, 5, 4],
             [5, 4, 4, ..., 4, 4, 5]])

    • lp
      (chain, draw)
      float64
      661.7 661.3 665.7 ... 665.6 668.8

      array([[661.70967519, 661.2875458 , 665.70729556, ..., 661.40453042,
              661.53660449, 661.51768248],
             [662.47186408, 665.25230302, 666.28391286, ..., 662.10849958,
              668.1822538 , 664.18500466],
             [662.60487008, 660.69369885, 662.99303071, ..., 663.00720061,
              664.60205279, 667.39059102],
             [664.1010232 , 660.27378818, 664.08879552, ..., 661.26338997,
              665.58839085, 668.77321037]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

  • Attributes: (4)

    created_at :

    2024-04-07T11:52:14.629957

    arviz_version :

    0.17.0

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:  (y_obs: 133)
  Coordinates:
    * y_obs    (y_obs) int64 0 1 2 3 4 5 6 7 8 ... 125 126 127 128 129 130 131 132
  Data variables:
      y        (y_obs) float64 0.0 -1.3 -2.7 0.0 -2.7 ... -2.7 10.7 -2.7 10.7
  Attributes:
      created_at:                  2024-04-07T11:52:14.631280
      arviz_version:               0.17.0
      inference_library:           numpyro
      inference_library_version:   0.13.2
      sampling_time:               1.871345
      modeling_interface:          bambi
      modeling_interface_version:  0.13.0

  xarray.Dataset

  • Dimensions:
    • y_obs: 133
  • Coordinates: (1)
    • y_obs
      (y_obs)
      int64
      0 1 2 3 4 5 ... 128 129 130 131 132

      array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
              14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
              28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
              42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
              56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
              70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
              84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
              98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
             112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
             126, 127, 128, 129, 130, 131, 132])

  • Data variables: (1)
    • y
      (y_obs)
      float64
      0.0 -1.3 -2.7 ... 10.7 -2.7 10.7

      array([   0. ,   -1.3,   -2.7,    0. ,   -2.7,   -2.7,   -2.7,   -1.3,
               -2.7,   -2.7,   -1.3,   -2.7,   -2.7,   -2.7,   -5.4,   -2.7,
               -5.4,    0. ,   -2.7,   -2.7,    0. ,  -13.3,   -5.4,   -5.4,
               -9.3,  -16. ,  -22.8,   -2.7,  -22.8,  -32.1,  -53.5,  -54.9,
              -40.2,  -21.5,  -21.5,  -50.8,  -42.9,  -26.8,  -21.5,  -50.8,
              -61.7,   -5.4,  -80.4,  -59. ,  -71. ,  -91.1,  -77.7,  -37.5,
              -85.6, -123.1, -101.9,  -99.1, -104.4, -112.5,  -50.8, -123.1,
              -85.6,  -72.3, -127.2, -123.1, -117.9, -134. , -101.9, -108.4,
             -123.1, -123.1, -128.5, -112.5,  -95.1,  -81.8,  -53.5,  -64.4,
              -57.6,  -72.3,  -44.3,  -26.8,   -5.4, -107.1,  -21.5,  -65.6,
              -16. ,  -45.6,  -24.2,    9.5,    4. ,   12. ,  -21.5,   37.5,
               46.9,  -17.4,   36.2,   75. ,    8.1,   54.9,   48.2,   46.9,
               16. ,   45.6,    1.3,   75. ,  -16. ,  -54.9,   69.6,   34.8,
               32.1,  -37.5,   22.8,   46.9,   10.7,    5.4,   -1.3,  -21.5,
              -13.3,   30.8,  -10.7,   29.4,    0. ,  -10.7,   14.7,   -1.3,
                0. ,   10.7,   10.7,  -26.8,  -14.7,  -13.3,    0. ,   10.7,
              -14.7,   -2.7,   10.7,   -2.7,   10.7])

  • Indexes: (1)
    • y_obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             123, 124, 125, 126, 127, 128, 129, 130, 131, 132],
            dtype='int64', name='y_obs', length=133))

  • Attributes: (7)

    created_at :

    2024-04-07T11:52:14.631280

    arviz_version :

    0.17.0

    inference_library :

    numpyro

    inference_library_version :

    0.13.2

    sampling_time :

    1.871345

    modeling_interface :

    bambi

    modeling_interface_version :

    0.13.0

  
  

and plot the model fit to see if it can recover the observed data.

ax = az.plot_hdi(new_data['X'], idata['posterior_predictive']['y'], fill_kwargs={'alpha': 0.2, 'color':'firebrick'}, figsize=(9, 8))

az.plot_hdi(new_data['X'], idata['posterior_predictive']['y'], fill_kwargs={'alpha': 0.8, 'color':'firebrick'}, hdi_prob=0.5)

y_mean = idata['posterior_predictive']['y'].mean(dim=('chain', 'draw'))

ax.plot(new_data['X'], y_mean, label='Expected posterior predictive', color='k')

ax.scatter(df['X'], df['y'], label='Observed Datapoints')

ax.legend()

ax.set_title("Posterior Predictive Distribution \n Based on HSGP approximation")

Text(0.5, 1.0, 'Posterior Predictive Distribution \n Based on HSGP approximation')

And we can compare versus the spline models to see that by the aggregate performance measures our HSGP model seems to come out on top.

models_dict = {"piecewise_constant": idata_spline1, "piecewise_linear": idata_spline2, "cubic_bspline": idata_spline3, "cubic_bspline_10": idata_spline4, 
"cubic_bspline_15": idata_spline5, 'hsgp': idata}
df_compare = az.compare(models_dict)
df_compare

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
hsgp	0	-608.819770	9.958827	0.000000	0.873213	12.009698	0.000000	False	log
cubic_bspline_10	1	-612.278059	11.465651	3.458289	0.000000	9.688433	2.962396	True	log
cubic_bspline_15	2	-616.252712	15.566839	7.432943	0.000000	9.833120	3.457579	True	log
cubic_bspline	3	-634.729956	8.698337	25.910186	0.000000	8.897139	7.686648	True	log
piecewise_constant	4	-643.781042	6.981098	34.961272	0.069506	9.770740	9.490817	False	log
piecewise_linear	5	-647.267410	6.170543	38.447641	0.057282	7.902853	9.200728	False	log

Recap

So far we’ve seen how we can use splines and gauassian processes to model highly eccentric functional relationships where the function could be approximated with univariate smoothing routine. Next we’ll show how to use hierarchical modelling over spline fits to extract insight into the data generating process over a family of curves. In particular we’ll focus on the development of insurance loss curves.

Insurance Loss Curves: Hierarchical Spline Models

We draw on car insurance losses data set discussed in Mick Cooney’s Stan case-study, but we simplify things for ourselves considerably by focusing on one type of loss and ensuring that each year under consideration has equal observations of the accruing losses.

loss_df = pd.read_csv('ppauto_pos.csv')
loss_df = loss_df[(loss_df['GRCODE'].isin([43, 353])) & (loss_df['DevelopmentYear'] < 1998)]

loss_df = loss_df[['GRCODE', 'AccidentYear', 'DevelopmentYear', 'DevelopmentLag', 'EarnedPremDIR_B', 'CumPaidLoss_B']]


loss_df.columns = ['grcode', 'acc_year', 'dev_year', 'dev_lag', 'premium', 'cum_loss']
loss_df['lr'] = loss_df['cum_loss'] / loss_df['premium']
loss_df = loss_df[(loss_df['acc_year'] <= 1992) & (loss_df['dev_lag'] <= 6)].reset_index(drop=True)

loss_df['year_code'] = loss_df['acc_year'].astype(str) + '_' + loss_df['grcode'].astype(str)
loss_df.sort_values(by=['year_code', 'acc_year', 'dev_lag'], inplace=True)
loss_df['standardised_premium'] = (loss_df.assign(mean_premium = np.mean(loss_df['premium']))
.assign(std_premium = np.std(loss_df['premium']))
.apply(lambda x: (x['mean_premium'] - x['premium']) /x['std_premium'], axis=1)
)

loss_df.head(12)

	grcode	acc_year	dev_year	dev_lag	premium	cum_loss	lr	year_code	standardised_premium
30	353	1988	1988	1	18793	4339	0.230884	1988_353	-0.347453
31	353	1988	1989	2	18793	9617	0.511733	1988_353	-0.347453
32	353	1988	1990	3	18793	11584	0.616400	1988_353	-0.347453
33	353	1988	1991	4	18793	12001	0.638589	1988_353	-0.347453
34	353	1988	1992	5	18793	12640	0.672591	1988_353	-0.347453
35	353	1988	1993	6	18793	12966	0.689938	1988_353	-0.347453
0	43	1988	1988	1	957	133	0.138976	1988_43	1.722344
1	43	1988	1989	2	957	333	0.347962	1988_43	1.722344
2	43	1988	1990	3	957	431	0.450366	1988_43	1.722344
3	43	1988	1991	4	957	570	0.595611	1988_43	1.722344
4	43	1988	1992	5	957	615	0.642633	1988_43	1.722344
5	43	1988	1993	6	957	615	0.642633	1988_43	1.722344

### Plot the Loss Curves Here we have plotted the developing loss curves from two different coded insurance products.

ax = (loss_df.pivot(index=['dev_lag'], columns=['grcode', 'acc_year'], values='lr')).plot(figsize=(10, 6))
ax.set_title("Loss Ratios by Year");

We want to model these curves collectively as instances of draws from a distribution of loss curves. To do so we will specify a PyMC hierarchical (mixed) spline model. To do so we will have a spline basis for the global hyper parameters beta_g and the individual parameters for each curve. Here we define a convenience function to generate the basis splines.

from patsy import bs as bs_patsy, dmatrix

def make_basis_splines(num_knots=3, max_dev=7):
    knot_list = np.linspace(0, max_dev, num_knots+2)[1:-1]
    dev_periods = np.arange(1, max_dev, 1)

    Bi = dmatrix(
        "bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True) - 1",
        {"dev_periods": dev_periods, "knots": knot_list},
    )

    Bg = dmatrix(
        "bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True) - 1",
        {"dev_periods": dev_periods, "knots": knot_list})


    return Bi, Bg

Bi, Bg = make_basis_splines()
Bg

DesignMatrix with shape (6, 7)
  Columns:
    ['bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[0]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[1]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[2]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[3]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[4]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[5]',
     'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)[6]']
  Terms:
    'bs_patsy(dev_periods, knots=knots, degree=3, include_intercept=True)' (columns 0:7)
  (to view full data, use np.asarray(this_obj))

Next we specify a model maker function to create the various pooled, unpooled and hierarhical (mixed) models of the insurance curve data. Note that even though we’re specifying a hierarhical model we have not specified a hierarchy over the insurance codes, instead we have added this as a “fixed” effect feature into our regression model.

def make_model(loss_df, num_knots=3, max_dev=7, model_type='mixed'):
    Bi, Bg = make_basis_splines(num_knots, max_dev)
    observed = loss_df['lr'].values
    uniques, unique_codes = pd.factorize(loss_df['year_code'])
    coords= {'years': unique_codes, 
            'splines': list(range(Bi.shape[1])) ,
            'measurement': list(range(6)), 
            'obs': uniques
            }

    with pm.Model(coords=coords) as sp_insur:
        basis_g = pm.MutableData('Bg', np.asfortranarray(Bg))

        tau = pm.HalfCauchy('tau', 1)
        ## Global Hierarchical Spline Terms
        beta_g = pm.Normal("beta_g", mu=0, sigma=tau, 
        dims='splines')
        mu_g = pm.Deterministic("mu_g", pm.math.dot(basis_g, beta_g), dims='measurement')

        ## Individual or Year Specific Spline Modifications
        if model_type in ['mixed', 'unpooled']:
            basis_i = pm.MutableData('Bi', np.asfortranarray(Bi))
            beta = pm.Normal("beta", mu=0, sigma=tau, dims=('splines', 'years'))
            mui = pm.Deterministic("mui", pm.math.dot(basis_i, beta), dims=('measurement', 'years'))
        
        ## Features
        prem = pm.MutableData('prem', loss_df['standardised_premium'].values)
        grcode = pm.MutableData('grcode', loss_df['grcode'] == 43)

        beta_prem = pm.Normal('beta_prem', 0, 1)
        beta_grcode = pm.Normal('beta_grcode', 0, 1)
        mu_prem = beta_prem*prem
        mu_grcode = beta_grcode*grcode

        ## Likelihood
        sigma = pm.TruncatedNormal("sigma", 1, lower=0.05)
        if model_type == 'mixed':
            mu = pm.Deterministic('mu',  mu_grcode + mu_prem + (mu_g.T + mui.T).ravel(), dims='obs')
            lr_likelihood = pm.Normal("lr", mu, sigma, observed=observed, dims=('obs'))
        elif model_type == 'pooled': 
             lr_likelihood = pm.Normal("lr",  np.repeat(mu_g, len(unique_codes)), sigma, observed=observed, dims='obs')
        elif model_type == 'unpooled':
            lr_likelihood = pm.Normal("lr",  mui.T.ravel(), sigma, observed=observed, dims=('obs'))


        ## Sampling
        idata_sp_insur = pm.sample(2000, return_inferencedata=True, target_accept=.99,
        idata_kwargs={"log_likelihood": True})
        idata_sp_insur = pm.sample_posterior_predictive(
            idata_sp_insur,extend_inferencedata=True)

    return idata_sp_insur, sp_insur


idata_sp_insur_unpooled, sp_insur_unpooled = make_model(loss_df, model_type='unpooled')
idata_sp_insur_pooled, sp_insur_pooled = make_model(loss_df, model_type='pooled')
idata_sp_insur_mixed, sp_insur_mixed = make_model(loss_df, model_type='mixed')

The model structure can be seen more clearly in this graph

pm.model_to_graphviz(sp_insur_mixed)

We can extract the effect of the differences grcodes and examine the baseline and annual spline related coefficients.

summary = az.summary(idata_sp_insur_mixed, var_names=['beta_grcode', 'beta_prem', 'beta_g', 'beta'])

summary

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
beta_grcode	0.330	0.074	0.195	0.469	0.001	0.001	4022.0	4343.0	1.0
beta_prem	-0.055	0.039	-0.131	0.016	0.001	0.000	4394.0	4646.0	1.0
beta_g[0]	0.075	0.088	-0.091	0.237	0.002	0.001	2979.0	4035.0	1.0
beta_g[1]	0.053	0.178	-0.281	0.387	0.002	0.002	7456.0	6406.0	1.0
beta_g[2]	0.439	0.136	0.193	0.697	0.002	0.001	5605.0	5857.0	1.0
...	...	...	...	...	...	...	...	...	...
beta[6, 1990_43]	0.550	0.103	0.343	0.733	0.002	0.001	3984.0	4846.0	1.0
beta[6, 1991_353]	-0.016	0.100	-0.211	0.166	0.002	0.001	3706.0	4830.0	1.0
beta[6, 1991_43]	0.107	0.103	-0.096	0.291	0.002	0.001	3992.0	4815.0	1.0
beta[6, 1992_353]	-0.030	0.100	-0.217	0.161	0.002	0.001	3992.0	5160.0	1.0
beta[6, 1992_43]	-0.067	0.126	-0.307	0.165	0.002	0.001	4186.0	5253.0	1.0

79 rows × 9 columns

Again we can compare the performance metrics of the various models.

compare_df = az.compare({'unpooled': idata_sp_insur_unpooled, 
            'pooled': idata_sp_insur_pooled, 
            'mixed': idata_sp_insur_mixed})

compare_df

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
mixed	0	63.898564	34.884147	0.000000	1.000000e+00	2.188274	0.000000	True	log
unpooled	1	42.518091	44.888145	21.380473	2.186085e-12	1.298397	1.538506	True	log
pooled	2	-9.214606	5.907311	73.113170	0.000000e+00	5.378668	4.429981	False	log

Plot the Posterior Predictive Checks

We can check how well the model can recapture the observed data.

def plot_ppc_splines(idata):
    fig, axs = plt.subplots(5, 2, figsize=(9, 30), sharey=True)
    axs = axs.flatten()
    dev_periods = np.arange(1, 7, 1)
    uniques, unique_codes = pd.factorize(loss_df['year_code'])
    for y, c in zip(unique_codes, range(10)):
        az.plot_hdi(dev_periods, idata['posterior_predictive']['lr'].sel(obs=c), color='firebrick', ax=axs[c], fill_kwargs={'alpha': 0.2}, hdi_prob=.89)
        az.plot_hdi(dev_periods, idata['posterior_predictive']['lr'].sel(obs=c), color='firebrick', ax=axs[c], hdi_prob=0.5)
        axs[c].scatter(dev_periods, loss_df[(loss_df['year_code'] == y)]['lr'], color='k', label='Actual Loss Ratio')
        axs[c].plot(dev_periods, loss_df[(loss_df['year_code'] == y)]['lr'], color='k')
        axs[c].set_title(f"Posterior Predictive Loss Ratio \n 89% and 50% HDI: {y}")
        axs[c].legend()
        axs[c].set_xlabel('Default Lag');

plot_ppc_splines(idata_sp_insur_mixed)

Plot the Hierarchical Components

In the following plot we show similarly how to recapture the observed data, but additionally we decompose the structure of the model in each case and extract baseline forecasts which would be our guide to future loss-ratio development curves in lieu of any other information.

mu_g = idata_sp_insur_mixed.posterior.stack(draws=("chain", "draw"))["mu_g"]

mu_i = idata_sp_insur_mixed.posterior.stack(draws=("chain", "draw"))["mui"]

mu  = idata_sp_insur_mixed.posterior.stack(draws=("chain", "draw"))['mu']

beta_grcode = idata_sp_insur_mixed.posterior.stack(draws=("chain", "draw"))['beta_grcode']

dev_periods = np.arange(1, 7, 1)
uniques, unique_codes = pd.factorize(loss_df['year_code'])

mosaic = """
         AB
         CD
         EF
         GH
         IJ
         KK
"""
fig, axs = plt.subplot_mosaic(mosaic, sharey=True, 
figsize=(9, 30))
axs = [axs[k] for k in axs.keys()] 

mu_g_mean = mu_g.mean(dim='draws')
for y, c in zip(unique_codes, range(10)):
    group_effect = 0
    if '43' in y: 
        group_effect = beta_grcode.mean().item()
    mu_i_mean = mu_i.sel(years=y).mean(dim='draws')
    axs[c].plot(dev_periods, group_effect + mu_g_mean.values + mu_i_mean.values, label='Combined + E(grp_effect)', color='purple', linewidth=3.5)
    axs[c].plot(dev_periods, group_effect + mu_g_mean.values, label='E(Hierarchical Baseline)', color='red', linestyle='--')
    axs[c].plot(dev_periods,  group_effect + mu_i_mean.values, label='E(Year Specific Adjustment term)', color='blue', linestyle='--')
    axs[c].scatter(dev_periods, loss_df[(loss_df['year_code'] == y)]['lr'], color='k', label='Actual Loss Ratio')
    axs[c].legend()
    az.plot_hdi(dev_periods,mu.sel(obs=c).T  , ax=axs[c], color='firebrick', fill_kwargs={'alpha': 0.2})
    az.plot_hdi(dev_periods, mu.sel(obs=c).T , ax=axs[c], color='firebrick', fill_kwargs={'alpha': 0.5}, hdi_prob=.50)
    axs[c].set_title(f"Components for Year {y}")

axs[10].plot(dev_periods, mu_g_mean.values, label='E(Hierarchical Baseline)', color='black')
axs[10].plot(dev_periods, mu_g_mean.values + group_effect, label='E(Hierarchical Baseline) + E(grp_effect)', color='black', linestyle='--')
az.plot_hdi(dev_periods, mu_g.T.values, color='slateblue', ax=axs[10], fill_kwargs={'alpha': 0.2})
az.plot_hdi(dev_periods, mu_g.T.values + group_effect, color='magenta', ax=axs[10], fill_kwargs={'alpha': 0.2})
az.plot_hdi(dev_periods, mu_g.T.values, color='slateblue', ax=axs[10], hdi_prob=.5)
az.plot_hdi(dev_periods, mu_g.T.values  + group_effect, color='magenta', ax=axs[10], hdi_prob=.5)
axs[10].set_title("Baseline Forecast Loss Ratio \n with Uncertainty Intervals")
axs[10].legend();

This plot is a bit clunky, because we’re mixing expectations and posterior distributions over the parameters. The point is just to highlight the “compositional” structure of our model. A better way to interrogate the implications of the model is to “push” forward different data through the posterior predictive distribution and derive a kind of ceteris paribus rule for accrual of losses. We d

with sp_insur_mixed: 
    pm.set_data({'grcode': np.ones(len(loss_df)), 
    })
    idata_43 = pm.sample_posterior_predictive(idata_sp_insur_mixed, var_names=['lr'], extend_inferencedata =True)

with sp_insur_mixed: 
    pm.set_data({'grcode': np.zeros(len(loss_df))})
    idata_353 = pm.sample_posterior_predictive(idata_sp_insur_mixed, var_names=['lr'], extend_inferencedata=True)

idata_353

Sampling: [lr]

100.00% [8000/8000 00:00<00:00]

Sampling: [lr]

100.00% [8000/8000 00:00<00:00]

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:      (chain: 4, draw: 2000, splines: 7, years: 10, measurement: 6,
                    obs: 60)
  Coordinates:
    * chain        (chain) int64 0 1 2 3
    * draw         (draw) int64 0 1 2 3 4 5 6 ... 1994 1995 1996 1997 1998 1999
    * splines      (splines) int64 0 1 2 3 4 5 6
    * years        (years) <U8 '1988_353' '1988_43' ... '1992_353' '1992_43'
    * measurement  (measurement) int64 0 1 2 3 4 5
    * obs          (obs) int64 0 0 0 0 0 0 1 1 1 1 1 1 ... 8 8 8 8 8 8 9 9 9 9 9 9
  Data variables:
      beta_g       (chain, draw, splines) float64 0.1471 -0.03058 ... 0.5207
      beta         (chain, draw, splines, years) float64 0.1228 ... -0.04697
      beta_prem    (chain, draw) float64 -0.05328 -0.08479 ... -0.02593 -0.04951
      beta_grcode  (chain, draw) float64 0.4148 0.3667 0.3818 ... 0.351 0.3281
      tau          (chain, draw) float64 0.2581 0.2447 0.2714 ... 0.2225 0.2312
      sigma        (chain, draw) float64 0.05043 0.05917 ... 0.0508 0.05102
      mu_g         (chain, draw, measurement) float64 0.1471 0.2455 ... 0.5207
      mui          (chain, draw, measurement, years) float64 0.1228 ... -0.04697
      mu           (chain, draw, obs) float64 0.2884 0.5118 ... 0.8584 0.8796
  Attributes:
      created_at:                 2024-04-07T11:53:53.245580
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3
      sampling_time:              39.531437158584595
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • splines: 7
    • years: 10
    • measurement: 6
    • obs: 60
  • Coordinates: (6)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999])

    • splines
      (splines)
      int64
      0 1 2 3 4 5 6

      array([0, 1, 2, 3, 4, 5, 6])

    • years
      (years)
      <U8
      '1988_353' '1988_43' ... '1992_43'

      array(['1988_353', '1988_43', '1989_353', '1989_43', '1990_353', '1990_43',
             '1991_353', '1991_43', '1992_353', '1992_43'], dtype='<U8')

    • measurement
      (measurement)
      int64
      0 1 2 3 4 5

      array([0, 1, 2, 3, 4, 5])

    • obs
      (obs)
      int64
      0 0 0 0 0 0 1 1 ... 8 8 9 9 9 9 9 9

      array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9])

  • Data variables: (9)
    • beta_g
      (chain, draw, splines)
      float64
      0.1471 -0.03058 ... 0.6265 0.5207

      array([[[ 0.14709482, -0.03057593,  0.36536821, ...,  0.56493346,
                0.38236167,  0.48785557],
              [ 0.05777286, -0.10195874,  0.43932537, ...,  0.75085442,
                0.32729495,  0.6014333 ],
              [-0.01735409,  0.3001635 ,  0.34119037, ...,  0.54193947,
                0.47805692,  0.55306759],
              ...,
              [-0.01266002, -0.20526287,  0.41746324, ...,  0.44991549,
                0.62809486,  0.50800348],
              [ 0.16210933,  0.10083014,  0.44979706, ...,  0.73709657,
                0.42675834,  0.55709393],
              [ 0.21241866, -0.00800354,  0.49882641, ...,  0.78231053,
                0.47610904,  0.73899835]],
      
             [[ 0.09886926,  0.11571252,  0.36536006, ...,  0.68783015,
                0.46738804,  0.61162468],
              [ 0.09117834,  0.04295766,  0.43724491, ...,  0.61708167,
                0.66425011,  0.58203739],
              [ 0.15108078,  0.31178546,  0.22084359, ...,  0.47529161,
                0.41188356,  0.64387646],
      ...
              [ 0.05331221, -0.39563518,  0.71260244, ...,  0.97217592,
                0.24803403,  0.56743485],
              [-0.01271875,  0.31190979,  0.33024387, ...,  0.44059866,
                0.33758343,  0.58697983],
              [ 0.04387177,  0.17895096,  0.47676331, ...,  0.54265466,
                0.58692165,  0.63674187]],
      
             [[ 0.06974003, -0.07411388,  0.51564521, ...,  0.694966  ,
                0.39137873,  0.48193329],
              [ 0.02693294,  0.04168749,  0.64152885, ...,  0.88033382,
                0.1140353 ,  0.41634181],
              [ 0.03853168,  0.03718208,  0.69019295, ...,  0.91836254,
                0.09589772,  0.4508061 ],
              ...,
              [-0.06683101,  0.19756622,  0.26813658, ...,  0.46359298,
                0.45872619,  0.53256906],
              [ 0.12607105,  0.2484163 ,  0.33416223, ...,  0.77120224,
                0.35186675,  0.5722983 ],
              [ 0.10991657,  0.08633124,  0.39636114, ...,  0.59471806,
                0.62648373,  0.52065511]]])

    • beta
      (chain, draw, splines, years)
      float64
      0.1228 -0.3016 ... 0.06587 -0.04697

      array([[[[ 1.22762129e-01, -3.01616146e-01,  1.03328507e-01, ...,
                -2.99612643e-01,  7.49347048e-02, -4.15317215e-01],
               [ 2.43803110e-01, -7.99853645e-02,  1.00467025e-02, ...,
                -5.43254755e-02, -1.29075688e-01, -4.60906636e-01],
               [ 2.90720611e-01, -2.91339632e-01,  3.76200583e-01, ...,
                 7.95590888e-02,  3.21321727e-01, -6.57439456e-02],
               ...,
               [ 3.66557601e-01, -2.27389241e-01,  6.10870835e-02, ...,
                -1.89668109e-01,  6.79335983e-02, -2.86258542e-02],
               [-1.32639078e-01, -1.58552911e-02,  2.24523853e-01, ...,
                 2.45358478e-01,  6.15860986e-03, -1.81421413e-02],
               [ 1.74464776e-01, -1.67943358e-01,  1.42439141e-01, ...,
                 1.29397424e-01,  4.67956869e-02, -2.43976450e-02]],
      
              [[ 1.75194381e-01, -1.19807005e-01,  1.32153439e-01, ...,
                -2.14812730e-01,  1.09640212e-01, -2.22434300e-01],
               [-3.61611172e-02, -2.19159377e-02,  2.97560312e-01, ...,
                 2.36196915e-02,  1.82200743e-01,  1.26924578e-01],
               [ 3.98931731e-01, -1.84929505e-01,  1.24334132e-01, ...,
                -5.84767973e-02,  4.35883309e-02, -3.17769413e-01],
      ...
               [ 5.06565185e-03, -3.34325814e-01,  2.40110601e-03, ...,
                -8.01798000e-02, -1.87972832e-01, -2.56137025e-01],
               [-6.40621365e-02, -8.79283433e-02, -7.43262798e-02, ...,
                 7.96750657e-02,  2.62457811e-01, -2.81438714e-02],
               [ 7.46187548e-02, -2.65938717e-01,  7.01132780e-02, ...,
                 1.03360769e-01, -1.05078831e-01, -5.71860731e-02]],
      
              [[ 7.32903661e-02, -2.20374824e-01,  1.38468692e-01, ...,
                -1.42030142e-01,  1.13265365e-01, -1.51380025e-01],
               [-9.14593004e-02,  1.57661663e-01,  3.08988599e-01, ...,
                -1.18975380e-01,  2.89423201e-01, -1.04312879e-01],
               [ 3.20066077e-01, -3.22129008e-01,  5.34809426e-02, ...,
                 9.78291771e-02,  9.66028750e-03, -2.55202372e-01],
               ...,
               [-3.55580383e-02, -2.40339616e-01, -1.59060524e-01, ...,
                 1.69944588e-01, -3.15972423e-02, -2.66585089e-01],
               [ 2.60770367e-01, -1.00729980e-01,  1.00356333e-01, ...,
                 6.06977373e-02, -3.14504966e-01,  1.50451241e-02],
               [ 1.98314078e-01, -1.17591360e-01,  3.80506957e-04, ...,
                 1.29902412e-01,  6.58737582e-02, -4.69687727e-02]]]])

    • beta_prem
      (chain, draw)
      float64
      -0.05328 -0.08479 ... -0.04951

      array([[-0.05327835, -0.08479205, -0.0497838 , ...,  0.0069594 ,
              -0.06376418, -0.04575061],
             [-0.04386236, -0.02639949, -0.04499374, ..., -0.03342609,
              -0.07446503, -0.07054487],
             [-0.05221887, -0.05901248, -0.10638247, ...,  0.0038883 ,
              -0.12298434, -0.13789127],
             [-0.04853679, -0.04860611, -0.03704374, ..., -0.11710353,
              -0.02593249, -0.0495121 ]])

    • beta_grcode
      (chain, draw)
      float64
      0.4148 0.3667 ... 0.351 0.3281

      array([[0.41483434, 0.36673737, 0.38175634, ..., 0.36934417, 0.34202595,
              0.22618224],
             [0.26570033, 0.24695166, 0.40395738, ..., 0.34040929, 0.25786865,
              0.43915123],
             [0.32312129, 0.31947762, 0.41299143, ..., 0.37462128, 0.36442257,
              0.35427986],
             [0.41647477, 0.39752618, 0.3792883 , ..., 0.46717039, 0.35104425,
              0.32806639]])

    • tau
      (chain, draw)
      float64
      0.2581 0.2447 ... 0.2225 0.2312

      array([[0.25814139, 0.24470884, 0.27136289, ..., 0.28101819, 0.28288784,
              0.21771619],
             [0.21919974, 0.25022668, 0.30999641, ..., 0.24978512, 0.26034655,
              0.24417764],
             [0.29731923, 0.25629949, 0.28035934, ..., 0.29610801, 0.29091214,
              0.24931157],
             [0.28583208, 0.26137499, 0.24486428, ..., 0.26322459, 0.22254574,
              0.2312261 ]])

    • sigma
      (chain, draw)
      float64
      0.05043 0.05917 ... 0.0508 0.05102

      array([[0.05043011, 0.05916892, 0.06425761, ..., 0.05333176, 0.05535915,
              0.06135004],
             [0.05357822, 0.05101462, 0.07464524, ..., 0.061423  , 0.06493909,
              0.05024091],
             [0.07680686, 0.05067121, 0.05663275, ..., 0.05136396, 0.05380478,
              0.05232874],
             [0.05156384, 0.05076381, 0.05468464, ..., 0.05590763, 0.05079953,
              0.0510178 ]])

    • mu_g
      (chain, draw, measurement)
      float64
      0.1471 0.2455 ... 0.6046 0.5207

      array([[[ 0.14709482,  0.24552231,  0.38498535,  0.45651708,
                0.48644369,  0.48785557],
              [ 0.05777286,  0.24912143,  0.35963165,  0.46883419,
                0.55875574,  0.6014333 ],
              [-0.01735409,  0.37233406,  0.53690373,  0.60614557,
                0.54140082,  0.55306759],
              ...,
              [-0.01266002,  0.24568294,  0.49718101,  0.51756979,
                0.5212902 ,  0.50800348],
              [ 0.16210933,  0.33337853,  0.43042292,  0.53033279,
                0.59734924,  0.55709393],
              [ 0.21241866,  0.33315404,  0.47512484,  0.57556518,
                0.64376597,  0.73899835]],
      
             [[ 0.09886926,  0.28715985,  0.38106435,  0.49304323,
                0.57907834,  0.61162468],
              [ 0.09117834,  0.31844337,  0.4576602 ,  0.52496821,
                0.61243674,  0.58203739],
              [ 0.15108078,  0.27200232,  0.33610536,  0.42243597,
                0.44760438,  0.64387646],
      ...
              [ 0.05331221,  0.30912012,  0.46372342,  0.55798223,
                0.65536774,  0.56743485],
              [-0.01271875,  0.32967537,  0.35913935,  0.39637844,
                0.40036066,  0.58697983],
              [ 0.04387177,  0.36904539,  0.41253961,  0.43926189,
                0.5324085 ,  0.63674187]],
      
             [[ 0.06974003,  0.31053797,  0.42603369,  0.4913884 ,
                0.5561625 ,  0.48193329],
              [ 0.02693294,  0.42437885,  0.51908433,  0.59986066,
                0.58252114,  0.41634181],
              [ 0.03853168,  0.4532827 ,  0.55246417,  0.62952408,
                0.60116973,  0.4508061 ],
              ...,
              [-0.06683101,  0.2787926 ,  0.41597981,  0.48440133,
                0.47054288,  0.53256906],
              [ 0.12607105,  0.3200615 ,  0.41612504,  0.56362392,
                0.60054526,  0.5722983 ],
              [ 0.10991657,  0.32510132,  0.50840252,  0.58098359,
                0.60459036,  0.52065511]]])

    • mui
      (chain, draw, measurement, years)
      float64
      0.1228 -0.3016 ... 0.06587 -0.04697

      array([[[[ 1.22762129e-01, -3.01616146e-01,  1.03328507e-01, ...,
                -2.99612643e-01,  7.49347048e-02, -4.15317215e-01],
               [ 2.47768186e-01, -2.24667376e-01,  2.28440858e-01, ...,
                 5.39123846e-02,  1.40799036e-01, -2.09344099e-01],
               [ 1.77621293e-01, -2.78254944e-01,  2.04971153e-01, ...,
                 1.23454996e-01,  1.25546117e-01, -1.57392734e-01],
               [ 1.98708556e-01, -2.56112083e-01,  1.01925498e-01, ...,
                 3.67454048e-02,  4.22598619e-02, -1.39905806e-01],
               [ 1.74740374e-01, -1.68507554e-01,  1.15404886e-01, ...,
                -6.73385054e-03,  3.84166402e-02, -5.16091016e-02],
               [ 1.74464776e-01, -1.67943358e-01,  1.42439141e-01, ...,
                 1.29397424e-01,  4.67956869e-02, -2.43976450e-02]],
      
              [[ 1.75194381e-01, -1.19807005e-01,  1.32153439e-01, ...,
                -2.14812730e-01,  1.09640212e-01, -2.22434300e-01],
               [ 2.19543544e-01, -1.37488373e-01,  1.81326643e-01, ...,
                -2.22792623e-02,  9.92125831e-02, -1.34140676e-01],
               [ 1.79223903e-01, -1.83645690e-01,  1.39052183e-01, ...,
                 1.82500926e-03,  8.50298630e-02, -1.17551635e-01],
               [ 5.50586116e-02, -1.32114202e-01,  1.25722338e-01, ...,
      ...
                 9.65206333e-02,  4.41982868e-02, -6.00265619e-02],
               [ 1.04240757e-01, -2.55349374e-01,  4.12326413e-02, ...,
                 7.87769466e-02, -4.30568812e-02, -8.78544594e-02],
               [ 1.89388173e-03, -2.44479053e-01, -1.63441372e-02, ...,
                 5.98968000e-03, -2.00611206e-02, -1.47204830e-01],
               [ 7.46187548e-02, -2.65938717e-01,  7.01132780e-02, ...,
                 1.03360769e-01, -1.05078831e-01, -5.71860731e-02]],
      
              [[ 7.32903661e-02, -2.20374824e-01,  1.38468692e-01, ...,
                -1.42030142e-01,  1.13265365e-01, -1.51380025e-01],
               [ 1.54303978e-01, -1.60313348e-01,  1.36364703e-01, ...,
                 6.54181114e-03,  8.71440563e-02, -1.82386905e-01],
               [ 1.28786040e-01, -2.53635508e-01,  5.77396278e-02, ...,
                -5.34519057e-04, -2.68008441e-02, -1.45442810e-01],
               [ 1.37521542e-02, -2.29206292e-01, -1.58014807e-02, ...,
                 2.56624004e-02, -4.75154853e-02, -1.47889537e-01],
               [ 6.07345095e-02, -1.93528770e-01, -4.88939800e-02, ...,
                 1.02998684e-01, -1.22675708e-01, -1.52030236e-01],
               [ 1.98314078e-01, -1.17591360e-01,  3.80506957e-04, ...,
                 1.29902412e-01,  6.58737582e-02, -4.69687727e-02]]]])

    • mu
      (chain, draw, obs)
      float64
      0.2884 0.5118 ... 0.8584 0.8796

      array([[[0.28836869, 0.51180224, 0.5811184 , ..., 0.81517292,
               0.93339624, 0.96201957],
              [0.26242853, 0.49812627, 0.56831685, ..., 0.86125482,
               0.8773677 , 0.95945981],
              [0.25626668, 0.49526178, 0.5073608 , ..., 0.8763235 ,
               0.95190741, 0.98022826],
              ...,
              [0.1115067 , 0.41263958, 0.68738084, ..., 0.90315159,
               0.85490935, 0.86564777],
              [0.13119547, 0.47147687, 0.62047344, ..., 0.84304468,
               0.95149832, 0.83330147],
              [0.28207232, 0.54288486, 0.55326974, ..., 0.87639481,
               0.88984391, 0.89896803]],
      
             [[0.14037368, 0.5129037 , 0.60073337, ..., 0.9115604 ,
               0.8881417 , 0.97808591],
              [0.22782439, 0.48684127, 0.56646073, ..., 0.87610374,
               1.00153701, 0.87293414],
              [0.18162724, 0.46963229, 0.48622689, ..., 0.84275107,
               0.98289781, 0.93886142],
      ...
              [0.20082738, 0.42643225, 0.60699477, ..., 0.86127559,
               0.84095966, 0.88260063],
              [0.15087005, 0.4870238 , 0.60383287, ..., 0.86466012,
               0.90485676, 0.96210462],
              [0.27027477, 0.48663867, 0.57454798, ..., 0.88434335,
               0.88414564, 0.97019646]],
      
             [[0.27283602, 0.45556152, 0.66937789, ..., 0.87633646,
               0.8852948 , 0.85235212],
              [0.19177964, 0.49100378, 0.58278649, ..., 0.83132431,
               0.86267181, 0.86751313],
              [0.26542108, 0.58263393, 0.62586374, ..., 0.88611531,
               0.95149841, 0.9255499 ],
              ...,
              [0.23437806, 0.47760062, 0.59005033, ..., 0.86292496,
               0.95994012, 0.96147194],
              [0.22273945, 0.57456583, 0.64510424, ..., 0.86756681,
               0.84513778, 0.90690958],
              [0.20041009, 0.49660845, 0.65439172, ..., 0.83896907,
               0.85843514, 0.87956135]]])

  • Indexes: (6)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999],
            dtype='int64', name='draw', length=2000))

    • splines
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6], dtype='int64', name='splines'))

    • years
      PandasIndex

      PandasIndex(Index(['1988_353', '1988_43', '1989_353', '1989_43', '1990_353', '1990_43',
             '1991_353', '1991_43', '1992_353', '1992_43'],
            dtype='object', name='years'))

    • measurement
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5], dtype='int64', name='measurement'))

    • obs
      PandasIndex

      PandasIndex(Index([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9],
            dtype='int64', name='obs'))

  • Attributes: (6)

    created_at :

    2024-04-07T11:53:53.245580

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    sampling_time :

    39.531437158584595

    tuning_steps :

    1000

  
  
• posterior_predictive

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 2000, obs: 60)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 ... 1993 1994 1995 1996 1997 1998 1999
    * obs      (obs) int64 0 0 0 0 0 0 1 1 1 1 1 1 2 ... 7 8 8 8 8 8 8 9 9 9 9 9 9
  Data variables:
      lr       (chain, draw, obs) float64 0.2116 0.4905 0.623 ... 0.823 0.8547
  Attributes:
      created_at:                 2024-04-07T11:53:55.312027
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • obs: 60
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999])

    • obs
      (obs)
      int64
      0 0 0 0 0 0 1 1 ... 8 8 9 9 9 9 9 9

      array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9])

  • Data variables: (1)
    • lr
      (chain, draw, obs)
      float64
      0.2116 0.4905 ... 0.823 0.8547

      array([[[0.21163338, 0.49046972, 0.62300337, ..., 0.75694077,
               0.93130723, 0.99109248],
              [0.27243884, 0.50601853, 0.52016259, ..., 0.89105769,
               0.83647018, 1.03984969],
              [0.20868832, 0.54714993, 0.50044958, ..., 0.87974445,
               1.03272604, 1.01649697],
              ...,
              [0.02456701, 0.37471305, 0.69343383, ..., 0.87838438,
               0.83809325, 0.7999993 ],
              [0.21190903, 0.45848566, 0.63509262, ..., 0.734679  ,
               0.96447683, 0.84938521],
              [0.22217771, 0.47399699, 0.53985622, ..., 0.80489967,
               0.82627374, 0.88208028]],
      
             [[0.15222622, 0.50072728, 0.58680085, ..., 0.92470976,
               0.91214001, 0.98791088],
              [0.17777105, 0.43588389, 0.56053856, ..., 0.89244992,
               1.08490795, 0.89902103],
              [0.12930594, 0.43790437, 0.33418056, ..., 0.90239426,
               0.91213598, 1.02251686],
      ...
              [0.14527766, 0.46356242, 0.60782167, ..., 0.8753737 ,
               0.7967867 , 0.90388222],
              [0.13115307, 0.53043118, 0.60719917, ..., 0.81805536,
               0.83861685, 1.02902821],
              [0.26431633, 0.53156251, 0.62786113, ..., 0.95079337,
               0.94697168, 0.89506736]],
      
             [[0.35042882, 0.34831396, 0.68425196, ..., 1.01296868,
               0.88344507, 0.88933241],
              [0.2216432 , 0.47058745, 0.54920618, ..., 0.86134596,
               0.85700789, 0.85610395],
              [0.27804952, 0.5611832 , 0.59016773, ..., 0.9510927 ,
               0.89595097, 0.95951637],
              ...,
              [0.28483144, 0.52634063, 0.6661701 , ..., 0.84203893,
               0.92800917, 0.94195816],
              [0.1306437 , 0.62010569, 0.68711115, ..., 0.8332255 ,
               0.87937284, 0.90651062],
              [0.07821064, 0.51200158, 0.62223873, ..., 0.92673155,
               0.82299107, 0.85471052]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999],
            dtype='int64', name='draw', length=2000))

    • obs
      PandasIndex

      PandasIndex(Index([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9],
            dtype='int64', name='obs'))

  • Attributes: (4)

    created_at :

    2024-04-07T11:53:55.312027

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

  
  
• log_likelihood

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 2000, obs: 60)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 ... 1993 1994 1995 1996 1997 1998 1999
    * obs      (obs) int64 0 0 0 0 0 0 1 1 1 1 1 1 2 ... 7 8 8 8 8 8 8 9 9 9 9 9 9
  Data variables:
      lr       (chain, draw, obs) float64 1.419 2.068 1.824 ... 1.98 1.662 1.67
  Attributes:
      created_at:                 2024-04-07T11:53:53.504735
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
    • obs: 60
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999])

    • obs
      (obs)
      int64
      0 0 0 0 0 0 1 1 ... 8 8 9 9 9 9 9 9

      array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9])

  • Data variables: (1)
    • lr
      (chain, draw, obs)
      float64
      1.419 2.068 1.824 ... 1.662 1.67

      array([[[ 1.41855205,  2.06822729,  1.82350233, ...,  1.69054312,
                1.89546228,  1.79058375],
              [ 1.76630683,  1.88197824,  1.57823057, ...,  1.90769601,
                1.80899707,  1.73327234],
              [ 1.74789748,  1.79306347,  0.38617686, ...,  1.7895876 ,
                1.54511257,  1.44903632],
              ...,
              [-0.49290401,  0.28609482,  1.1265904 , ...,  1.66964594,
                1.5929081 ,  1.40465211],
              [ 0.35361235,  1.71057732,  1.97226719, ...,  1.93342651,
                1.60304303,  0.61979427],
              [ 1.52413635,  1.74330524,  1.3427872 , ...,  1.83203795,
                1.84652275,  1.78602787]],
      
             [[ 0.58079207,  2.00743538,  1.96492473, ...,  1.52654105,
                1.9652278 ,  1.50641378],
              [ 2.05490621,  1.93766417,  1.57756694, ...,  2.00051921,
                0.21965248,  1.54702466],
              [ 1.45835184,  1.51701527,  0.15549988, ...,  1.65236915,
                1.11396567,  1.65740746],
      ...
              [ 1.87867014,  0.67089452,  2.03311643, ...,  2.048901  ,
                1.30261907,  1.71811733],
              [ 0.89770171,  1.89800392,  1.97617846, ...,  1.99792668,
                2.00324389,  1.75844059],
              [ 1.74794795,  1.91628537,  1.71144291, ...,  1.91401833,
                1.96107665,  1.64898321]],
      
             [[ 1.71502615,  1.45264382,  1.51819123, ...,  1.98949435,
                1.98192954,  1.06874428],
              [ 1.76493928,  1.9782589 ,  1.84241207, ...,  1.91298981,
                1.73419019,  1.43285229],
              [ 1.78779273,  1.14672571,  1.972258  , ...,  1.86432585,
                1.60606894,  1.98702791],
              ...,
              [ 1.96316284,  1.77875141,  1.854053  , ...,  1.9626551 ,
                1.46009392,  1.74574659],
              [ 2.04807766,  1.29599819,  1.90128613, ...,  2.04671906,
                1.39525368,  2.00138456],
              [ 1.87824895,  2.01269842,  1.77936715, ...,  1.97954322,
                1.66213543,  1.66973209]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999],
            dtype='int64', name='draw', length=2000))

    • obs
      PandasIndex

      PandasIndex(Index([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9],
            dtype='int64', name='obs'))

  • Attributes: (4)

    created_at :

    2024-04-07T11:53:53.504735

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:                (chain: 4, draw: 2000)
  Coordinates:
    * chain                  (chain) int64 0 1 2 3
    * draw                   (draw) int64 0 1 2 3 4 5 ... 1995 1996 1997 1998 1999
  Data variables: (12/17)
      lp                     (chain, draw) float64 82.88 82.71 ... 89.52 83.05
      index_in_trajectory    (chain, draw) int64 -103 -79 -103 59 ... -67 60 84 50
      acceptance_rate        (chain, draw) float64 0.9892 0.972 ... 0.9727 0.9946
      energy                 (chain, draw) float64 -45.6 -45.6 ... -51.22 -49.21
      step_size              (chain, draw) float64 0.03994 0.03994 ... 0.03634
      n_steps                (chain, draw) float64 127.0 127.0 ... 127.0 127.0
      ...                     ...
      tree_depth             (chain, draw) int64 7 7 7 7 7 7 7 8 ... 7 8 6 7 7 7 7
      perf_counter_start     (chain, draw) float64 1.549e+06 ... 1.549e+06
      perf_counter_diff      (chain, draw) float64 0.009798 0.009783 ... 0.009734
      energy_error           (chain, draw) float64 -0.01187 0.0464 ... -0.04328
      diverging              (chain, draw) bool False False False ... False False
      smallest_eigval        (chain, draw) float64 nan nan nan nan ... nan nan nan
  Attributes:
      created_at:                 2024-04-07T11:53:53.257993
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3
      sampling_time:              39.531437158584595
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 2000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 ... 1996 1997 1998 1999

      array([   0,    1,    2, ..., 1997, 1998, 1999])

  • Data variables: (17)
    • lp
      (chain, draw)
      float64
      82.88 82.71 70.44 ... 89.52 83.05

      array([[82.88314656, 82.70820473, 70.43635026, ..., 64.04408313,
              75.61252313, 76.60720398],
             [75.23468291, 75.03510548, 52.31977537, ..., 70.0051281 ,
              74.63419241, 83.31401036],
             [62.58771038, 73.64498621, 68.42080124, ..., 69.88669862,
              80.27394702, 88.14011222],
             [77.76022401, 81.46355219, 85.76862643, ..., 83.6520337 ,
              89.52369819, 83.04884813]])

    • index_in_trajectory
      (chain, draw)
      int64
      -103 -79 -103 59 ... -67 60 84 50

      array([[-103,  -79, -103, ...,   84, -107,  -61],
             [  48,   35,   86, ...,   76, -217,   84],
             [  76,   38,   58, ...,   54, -112,  -49],
             [  68,  -44,  -23, ...,   60,   84,   50]])

    • acceptance_rate
      (chain, draw)
      float64
      0.9892 0.972 ... 0.9727 0.9946

      array([[0.98915436, 0.97202427, 0.98573292, ..., 0.99522913, 0.99918183,
              0.98222681],
             [0.9572357 , 0.93792687, 0.972995  , ..., 0.9838465 , 0.99659148,
              0.99459224],
             [0.94662177, 0.99894729, 0.99917357, ..., 0.99231341, 0.99976952,
              0.97727941],
             [0.99774905, 0.99742566, 0.99997425, ..., 0.9761266 , 0.97267422,
              0.9946227 ]])

    • energy
      (chain, draw)
      float64
      -45.6 -45.6 ... -51.22 -49.21

      array([[-45.59915423, -45.60028451, -39.27953877, ..., -21.1341642 ,
              -26.23847221, -36.71135969],
             [-31.89534318, -37.5928412 , -19.65206341, ..., -39.33146564,
              -26.85954691, -44.17277685],
             [-15.93298176, -31.38144429, -37.930313  , ..., -31.10975541,
              -39.49420063, -46.77832581],
             [-35.12324858, -46.26656322, -55.15999839, ..., -47.68734374,
              -51.2180004 , -49.21403014]])

    • step_size
      (chain, draw)
      float64
      0.03994 0.03994 ... 0.03634 0.03634

      array([[0.03994314, 0.03994314, 0.03994314, ..., 0.03994314, 0.03994314,
              0.03994314],
             [0.03539491, 0.03539491, 0.03539491, ..., 0.03539491, 0.03539491,
              0.03539491],
             [0.03785024, 0.03785024, 0.03785024, ..., 0.03785024, 0.03785024,
              0.03785024],
             [0.03633869, 0.03633869, 0.03633869, ..., 0.03633869, 0.03633869,
              0.03633869]])

    • n_steps
      (chain, draw)
      float64
      127.0 127.0 127.0 ... 127.0 127.0

      array([[127., 127., 127., ..., 127., 255., 255.],
             [255., 127., 127., ..., 127., 639., 127.],
             [127., 127., 127., ..., 127., 127., 127.],
             [127., 127., 127., ..., 127., 127., 127.]])

    • process_time_diff
      (chain, draw)
      float64
      0.009776 0.009782 ... 0.009725

      array([[0.009776, 0.009782, 0.009781, ..., 0.009794, 0.01969 , 0.019518],
             [0.019214, 0.009486, 0.009483, ..., 0.009694, 0.048177, 0.009784],
             [0.009565, 0.009767, 0.009734, ..., 0.009696, 0.009687, 0.009718],
             [0.009759, 0.009842, 0.009825, ..., 0.009774, 0.009824, 0.009725]])

    • max_energy_error
      (chain, draw)
      float64
      0.0478 0.0785 ... 0.07684 -0.05871

      array([[ 0.04780167,  0.07850324,  0.06061101, ..., -0.05307146,
              -0.04432925,  0.06313809],
             [ 0.08892657,  0.15311003,  0.10026589, ...,  0.0389184 ,
              -0.04174343, -0.05848511],
             [ 0.16877403, -0.09043713, -0.05569405, ...,  0.02686471,
              -0.02740915,  0.05242059],
             [-0.07282855, -0.03674812, -0.05754739, ...,  0.06201199,
               0.07684371, -0.05871458]])

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • step_size_bar
      (chain, draw)
      float64
      0.03811 0.03811 ... 0.04005 0.04005

      array([[0.03810643, 0.03810643, 0.03810643, ..., 0.03810643, 0.03810643,
              0.03810643],
             [0.03717455, 0.03717455, 0.03717455, ..., 0.03717455, 0.03717455,
              0.03717455],
             [0.03980489, 0.03980489, 0.03980489, ..., 0.03980489, 0.03980489,
              0.03980489],
             [0.04004681, 0.04004681, 0.04004681, ..., 0.04004681, 0.04004681,
              0.04004681]])

    • tree_depth
      (chain, draw)
      int64
      7 7 7 7 7 7 7 8 ... 7 7 8 6 7 7 7 7

      array([[ 7,  7,  7, ...,  7,  8,  8],
             [ 8,  7,  7, ...,  7, 10,  7],
             [ 7,  7,  7, ...,  7,  7,  7],
             [ 7,  7,  7, ...,  7,  7,  7]])

    • perf_counter_start
      (chain, draw)
      float64
      1.549e+06 1.549e+06 ... 1.549e+06

      array([[1548942.31224196, 1548942.32212742, 1548942.33198654, ...,
              1548967.86465525, 1548967.87452525, 1548967.8943025 ],
             [1548941.62282217, 1548941.64214321, 1548941.65169504, ...,
              1548967.154324  , 1548967.16408525, 1548967.21236642],
             [1548941.26976637, 1548941.27945367, 1548941.28929225, ...,
              1548962.55801167, 1548962.56777346, 1548962.57752988],
             [1548941.39335158, 1548941.40319771, 1548941.41314825, ...,
              1548963.33556242, 1548963.34542987, 1548963.35536142]])

    • perf_counter_diff
      (chain, draw)
      float64
      0.009798 0.009783 ... 0.009734

      array([[0.00979821, 0.00978275, 0.00978029, ..., 0.0097945 , 0.01969233,
              0.01951783],
             [0.01921813, 0.00948567, 0.00948379, ..., 0.00969371, 0.04819417,
              0.009784  ],
             [0.00957462, 0.00976708, 0.00973396, ..., 0.00969692, 0.00968704,
              0.00973129],
             [0.00975842, 0.00985158, 0.009826  , ..., 0.00978396, 0.00985   ,
              0.00973404]])

    • energy_error
      (chain, draw)
      float64
      -0.01187 0.0464 ... 0.0433 -0.04328

      array([[-0.01187417,  0.04640385,  0.04550437, ..., -0.00464313,
              -0.02107491,  0.03507634],
             [ 0.04082534,  0.0533617 ,  0.05805446, ...,  0.02935723,
              -0.00539725,  0.02254652],
             [ 0.02193904, -0.05320018,  0.01166856, ..., -0.00132098,
              -0.004021  ,  0.03337353],
             [-0.06111826,  0.00941168, -0.02694765, ...,  0.00924557,
               0.04330169, -0.04327881]])

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
             ...
             1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999],
            dtype='int64', name='draw', length=2000))

  • Attributes: (6)

    created_at :

    2024-04-07T11:53:53.257993

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

    sampling_time :

    39.531437158584595

    tuning_steps :

    1000

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:  (obs: 60)
  Coordinates:
    * obs      (obs) int64 0 0 0 0 0 0 1 1 1 1 1 1 2 ... 7 8 8 8 8 8 8 9 9 9 9 9 9
  Data variables:
      lr       (obs) float64 0.2309 0.5117 0.6164 0.6386 ... 0.859 0.9038 0.9244
  Attributes:
      created_at:                 2024-04-07T11:53:53.260548
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3

  xarray.Dataset

  • Dimensions:
    • obs: 60
  • Coordinates: (1)
    • obs
      (obs)
      int64
      0 0 0 0 0 0 1 1 ... 8 8 9 9 9 9 9 9

      array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9])

  • Data variables: (1)
    • lr
      (obs)
      float64
      0.2309 0.5117 ... 0.9038 0.9244

      array([0.23088384, 0.51173309, 0.61639972, 0.63858884, 0.67259086,
             0.68993774, 0.13897597, 0.34796238, 0.45036573, 0.59561129,
             0.64263323, 0.64263323, 0.2345894 , 0.47973401, 0.58713321,
             0.62354866, 0.63473718, 0.63742875, 0.25277402, 0.47253045,
             0.64005413, 0.69797023, 0.74776725, 0.80270636, 0.29829006,
             0.51463926, 0.59872363, 0.66151898, 0.68509768, 0.70209967,
             0.33072662, 0.79244053, 1.12088628, 1.31753014, 1.40029326,
             1.42440534, 0.25679105, 0.44820942, 0.53957139, 0.57970674,
             0.59681361, 0.59784754, 0.25876918, 0.65744596, 0.86254492,
             0.94998004, 0.98790852, 1.03097017, 0.24606574, 0.41036624,
             0.49781029, 0.54771644, 0.56508975, 0.57923865, 0.25779626,
             0.54739102, 0.76565216, 0.85900276, 0.90375243, 0.9244402 ])

  • Indexes: (1)
    • obs
      PandasIndex

      PandasIndex(Index([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
             4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7,
             8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9],
            dtype='int64', name='obs'))

  • Attributes: (4)

    created_at :

    2024-04-07T11:53:53.260548

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

  
  
• constant_data

  <xarray.Dataset>
  Dimensions:       (Bg_dim_0: 6, Bg_dim_1: 7, Bi_dim_0: 6, Bi_dim_1: 7,
                     prem_dim_0: 60, grcode_dim_0: 60)
  Coordinates:
    * Bg_dim_0      (Bg_dim_0) int64 0 1 2 3 4 5
    * Bg_dim_1      (Bg_dim_1) int64 0 1 2 3 4 5 6
    * Bi_dim_0      (Bi_dim_0) int64 0 1 2 3 4 5
    * Bi_dim_1      (Bi_dim_1) int64 0 1 2 3 4 5 6
    * prem_dim_0    (prem_dim_0) int64 0 1 2 3 4 5 6 7 ... 52 53 54 55 56 57 58 59
    * grcode_dim_0  (grcode_dim_0) int64 0 1 2 3 4 5 6 7 ... 53 54 55 56 57 58 59
  Data variables:
      Bg            (Bg_dim_0, Bg_dim_1) float64 1.0 0.0 0.0 0.0 ... 0.0 0.0 1.0
      Bi            (Bi_dim_0, Bi_dim_1) float64 1.0 0.0 0.0 0.0 ... 0.0 0.0 1.0
      prem          (prem_dim_0) float64 -0.3475 -0.3475 -0.3475 ... -1.572 -1.572
      grcode        (grcode_dim_0) float64 0.0 0.0 0.0 0.0 0.0 ... 1.0 1.0 1.0 1.0
  Attributes:
      created_at:                 2024-04-07T11:53:53.261451
      arviz_version:              0.17.0
      inference_library:          pymc
      inference_library_version:  5.10.3

  xarray.Dataset

  • Dimensions:
    • Bg_dim_0: 6
    • Bg_dim_1: 7
    • Bi_dim_0: 6
    • Bi_dim_1: 7
    • prem_dim_0: 60
    • grcode_dim_0: 60
  • Coordinates: (6)
    • Bg_dim_0
      (Bg_dim_0)
      int64
      0 1 2 3 4 5

      array([0, 1, 2, 3, 4, 5])

    • Bg_dim_1
      (Bg_dim_1)
      int64
      0 1 2 3 4 5 6

      array([0, 1, 2, 3, 4, 5, 6])

    • Bi_dim_0
      (Bi_dim_0)
      int64
      0 1 2 3 4 5

      array([0, 1, 2, 3, 4, 5])

    • Bi_dim_1
      (Bi_dim_1)
      int64
      0 1 2 3 4 5 6

      array([0, 1, 2, 3, 4, 5, 6])

    • prem_dim_0
      (prem_dim_0)
      int64
      0 1 2 3 4 5 6 ... 54 55 56 57 58 59

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59])

    • grcode_dim_0
      (grcode_dim_0)
      int64
      0 1 2 3 4 5 6 ... 54 55 56 57 58 59

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59])

  • Data variables: (4)
    • Bg
      (Bg_dim_0, Bg_dim_1)
      float64
      1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 1.0

      array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 3.08571429e-01, 5.69339736e-01, 1.21488595e-01,
              6.00240096e-04, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 1.14285714e-02, 4.09819928e-01, 5.03721489e-01,
              7.50300120e-02, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 7.50300120e-02, 5.03721489e-01,
              4.09819928e-01, 1.14285714e-02, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 6.00240096e-04, 1.21488595e-01,
              5.69339736e-01, 3.08571429e-01, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])

    • Bi
      (Bi_dim_0, Bi_dim_1)
      float64
      1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 1.0

      array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 3.08571429e-01, 5.69339736e-01, 1.21488595e-01,
              6.00240096e-04, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 1.14285714e-02, 4.09819928e-01, 5.03721489e-01,
              7.50300120e-02, 0.00000000e+00, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 7.50300120e-02, 5.03721489e-01,
              4.09819928e-01, 1.14285714e-02, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 6.00240096e-04, 1.21488595e-01,
              5.69339736e-01, 3.08571429e-01, 0.00000000e+00],
             [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
              0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])

    • prem
      (prem_dim_0)
      float64
      -0.3475 -0.3475 ... -1.572 -1.572

      array([-0.34745346, -0.34745346, -0.34745346, -0.34745346, -0.34745346,
             -0.34745346,  1.72234378,  1.72234378,  1.72234378,  1.72234378,
              1.72234378,  1.72234378, -0.3654406 , -0.3654406 , -0.3654406 ,
             -0.3654406 , -0.3654406 , -0.3654406 ,  1.40460971,  1.40460971,
              1.40460971,  1.40460971,  1.40460971,  1.40460971, -0.5486773 ,
             -0.5486773 , -0.5486773 , -0.5486773 , -0.5486773 , -0.5486773 ,
              1.12110923,  1.12110923,  1.12110923,  1.12110923,  1.12110923,
              1.12110923, -0.63582788, -0.63582788, -0.63582788, -0.63582788,
             -0.63582788, -0.63582788, -0.20123544, -0.20123544, -0.20123544,
             -0.20123544, -0.20123544, -0.20123544, -0.5779209 , -0.5779209 ,
             -0.5779209 , -0.5779209 , -0.5779209 , -0.5779209 , -1.57150713,
             -1.57150713, -1.57150713, -1.57150713, -1.57150713, -1.57150713])

    • grcode
      (grcode_dim_0)
      float64
      0.0 0.0 0.0 0.0 ... 1.0 1.0 1.0 1.0

      array([0., 0., 0., 0., 0., 0., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0.,
             0., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1., 1.,
             1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1., 1., 1., 1., 0., 0., 0.,
             0., 0., 0., 1., 1., 1., 1., 1., 1.])

  • Indexes: (6)
    • Bg_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5], dtype='int64', name='Bg_dim_0'))

    • Bg_dim_1
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6], dtype='int64', name='Bg_dim_1'))

    • Bi_dim_0
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5], dtype='int64', name='Bi_dim_0'))

    • Bi_dim_1
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3, 4, 5, 6], dtype='int64', name='Bi_dim_1'))

    • prem_dim_0
      PandasIndex

      PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59],
            dtype='int64', name='prem_dim_0'))

    • grcode_dim_0
      PandasIndex

      PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59],
            dtype='int64', name='grcode_dim_0'))

  • Attributes: (4)

    created_at :

    2024-04-07T11:53:53.261451

    arviz_version :

    0.17.0

    inference_library :

    pymc

    inference_library_version :

    5.10.3

  
  

Even here though we want to average the curves over the specific years in the data and abstract a view of the model implications under different counterfactual settings. Here we define a helper function to effect this step.

def get_posterior_predictive_curve(idata, prem=2, grcode=1):
    weighted_splines = [np.dot(np.asfortranarray(Bi), az.extract(idata['posterior']['beta'])['beta'].values[:, :, i]) for i in range(1000)]

    weighted_splines_1 = [np.dot(np.asfortranarray(Bg), az.extract(idata['posterior']['beta_g'])['beta_g'].values[:, i]) for i in range(1000)]

    beta_grcode = az.extract(idata['posterior']['beta_grcode'])['beta_grcode']

    beta_prem = az.extract(idata['posterior']['beta_prem'])['beta_prem']
    df1 = pd.DataFrame([beta_prem.values[i]*prem + beta_grcode.values[i]*grcode for i in range(1000)]).T

    ## Crucial step averaging over the years to get
    ## a view of the development period
    df = pd.concat([pd.DataFrame(weighted_splines_1[i].T + weighted_splines[i].T).mean() for i in range(1000)], axis=1)

    df = df1.iloc[0] + df

    return df

pred_df_1 = get_posterior_predictive_curve(idata_43, prem=1, grcode=1)

pred_df_0 = get_posterior_predictive_curve(idata_353, prem=1, grcode=0) 

fig, ax = plt.subplots(figsize=(9, 7), sharey=True)


ax.plot(pred_df_1, color='firebrick', alpha=0.1)
ax.plot(pred_df_0, color='slateblue', alpha=0.2);
ax.plot(pred_df_0.mean(axis=1), linestyle='-', color='k', linewidth=4, label='grcode 353 prem 1')

ax.plot(pred_df_1.mean(axis=1), linestyle='--', color='grey', linewidth=4, label='grcode 43 prem 1')

ax.set_title("Posterior Samples of the Trajectories \n Under different Counterfactual settings")
ax.set_ylabel("Loss Ratio")
ax.set_xlabel("Development Period")
ax.legend();

Multiple Smoother Regression Models

To demonstrate how spline modelling can be further adapting to the multiple regression like cases we use the PISA data set discussed in this case study from Michael Clark.

The data set has been constructed using average Science scores by country from the Programme for International Student Assessment (PISA) 2006, along with GNI per capita (Purchasing Power Parity, 2005 dollars), Educational Index, Health Index, and Human Development Index from UN data. We want to model the overall outcome score as a function of these broad demographic features.

pisa_df = pd.read_csv("https://raw.githubusercontent.com/m-clark/generalized-additive-models/master/data/pisasci2006.csv")

pisa_df

	Country	Overall	Issues	Explain	Evidence	Interest	Support	Income	Health	Edu	HDI
0	Albania	NaN	NaN	NaN	NaN	NaN	NaN	0.599	0.886	0.716	0.724
1	Argentina	391.0	395.0	386.0	385.0	567.0	506.0	0.678	0.868	0.786	0.773
2	Australia	527.0	535.0	520.0	531.0	465.0	487.0	0.826	0.965	0.978	0.920
3	Austria	511.0	505.0	516.0	505.0	507.0	515.0	0.835	0.944	0.824	0.866
4	Azerbaijan	382.0	353.0	412.0	344.0	612.0	542.0	0.566	0.780	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...
60	Turkey	424.0	427.0	423.0	417.0	540.0	563.0	0.679	0.828	0.562	0.681
61	United Arab Emirates	NaN	NaN	NaN	NaN	NaN	NaN	0.909	0.878	0.686	0.818
62	United Kingdom	515.0	514.0	517.0	514.0	464.0	470.0	0.833	0.934	0.798	0.853
63	United States	489.0	492.0	486.0	489.0	480.0	490.0	0.872	0.911	0.930	0.904
64	Uruguay	428.0	429.0	423.0	429.0	567.0	510.0	0.658	0.884	0.740	0.755

65 rows × 11 columns

The relationships displayed between each of these measures is not obviously linear, and as such could plausibly benefit from being modelled with splines.

g = sns.pairplot(pisa_df[['Overall', 'Income', 'Support', 'Health', 'Edu']],  kind="reg")
g.fig.suptitle("Pair Plot of Complex Relations", y=1.05);

We define two basic models for contrasting. Note here how we have to define a seperate spline basis for each of the covariates.

formula = "Overall ~ Income + Edu + Health"
base_model = bmb.Model(formula, pisa_df, dropna=True)

knots_income = np.linspace(np.min(pisa_df['Income']), np.max(pisa_df['Income']), 5+2)[1:-1]

knots_edu = np.linspace(np.min(pisa_df['Edu']), np.max(pisa_df['Edu']), 5+2)[1:-1]

knots_health = np.linspace(np.min(pisa_df['Health']), np.max(pisa_df['Health']), 5+2)[1:-1]

formula_spline = """Overall ~ bs(Income, degree=3, knots=knots_income) + bs(Edu, degree=3, knots=knots_edu) + bs(Health, degree=3, knots=knots_health)"""

spline_model = bmb.Model(formula_spline, pisa_df, dropna=True)

spline_model

Automatically removing 13/65 rows from the dataset.

Automatically removing 13/65 rows from the dataset.

       Formula: Overall ~ bs(Income, degree=3, knots=knots_income) + bs(Edu, degree=3, knots=knots_edu) + bs(Health, degree=3, knots=knots_health)
        Family: gaussian
          Link: mu = identity
  Observations: 52
        Priors: 
    target = mu
        Common-level effects
            Intercept ~ Normal(mu: 471.1538, sigma: 464.5641)
            bs(Income, degree=3, knots=knots_income) ~ Normal(mu: [0. 0. 0. 0. 0. 0. 0. 0.], sigma:
                [2700.1849 1474.917   742.9386  517.2898  651.4592  515.7397  867.8945
              876.0978])
            bs(Edu, degree=3, knots=knots_edu) ~ Normal(mu: [0. 0. 0. 0. 0. 0. 0. 0.], sigma: [1120.7928
                920.8213  702.8955  537.103   560.453   582.8883 1023.4867
              680.9089])
            bs(Health, degree=3, knots=knots_health) ~ Normal(mu: [0. 0. 0. 0. 0. 0. 0. 0.], sigma:
                [1465.0684 3939.0308  690.5161  531.6981  655.7589  552.7209  642.9461
              925.9747])
        
        Auxiliary parameters
            sigma ~ HalfStudentT(nu: 4.0, sigma: 53.4654)

spline_model.build()
spline_model.graph()

base_idata = base_model.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True})
spline_idata = spline_model.fit(random_seed=random_seed, idata_kwargs={"log_likelihood": True}, target_accept=.95)
base_model.predict(base_idata, kind='pps')
spline_model.predict(spline_idata, kind='pps')

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [Overall_sigma, Intercept, Income, Edu, Health]

100.00% [8000/8000 00:02<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [Overall_sigma, Intercept, bs(Income, degree=3, knots=knots_income), bs(Edu, degree=3, knots=knots_edu), bs(Health, degree=3, knots=knots_health)]

100.00% [8000/8000 00:27<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 27 seconds.

We can compare the simple regression approach to the spline based regression in the now usual way.

compare_df = az.compare({'spline': spline_idata, 'raw': base_idata})

compare_df

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
spline	0	-254.195127	20.312286	0.000000	0.603434	4.737638	0.000000	True	log
raw	1	-263.049613	9.606047	8.854486	0.396566	10.073402	9.899356	True	log

The coefficients comparisons are harder

az.summary(base_idata)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
Intercept	118.438	81.996	-41.166	266.033	1.522	1.086	2915.0	2986.0	1.0
Income	181.955	89.547	17.413	353.598	1.913	1.353	2196.0	2432.0	1.0
Edu	234.125	59.153	120.739	340.343	1.017	0.740	3414.0	2518.0	1.0
Health	30.395	140.966	-223.699	298.258	2.852	2.153	2450.0	2598.0	1.0
Overall_sigma	34.234	3.620	27.882	41.175	0.060	0.043	3757.0	2879.0	1.0
Overall_mean[0]	452.208	6.052	440.799	463.946	0.092	0.065	4309.0	2865.0	1.0
Overall_mean[1]	527.038	10.136	508.570	546.504	0.176	0.125	3315.0	3356.0	1.0
Overall_mean[2]	491.982	6.878	479.124	505.280	0.117	0.083	3444.0	3314.0	1.0
Overall_mean[3]	501.282	6.069	490.053	512.985	0.101	0.072	3600.0	3132.0	1.0
Overall_mean[4]	410.451	8.866	394.136	427.611	0.135	0.095	4339.0	3107.0	1.0
Overall_mean[5]	446.420	6.947	433.537	459.774	0.101	0.071	4778.0	3353.0	1.0
Overall_mean[6]	511.366	6.848	498.881	525.124	0.111	0.078	3826.0	3258.0	1.0
Overall_mean[7]	447.819	11.282	426.627	468.792	0.221	0.156	2601.0	2310.0	1.0
Overall_mean[8]	481.644	10.339	460.677	499.738	0.173	0.123	3571.0	2874.0	1.0
Overall_mean[9]	401.813	9.440	383.890	419.274	0.141	0.099	4526.0	3355.0	1.0
Overall_mean[10]	455.263	5.193	446.072	465.558	0.076	0.054	4743.0	2812.0	1.0
Overall_mean[11]	501.446	8.079	486.819	517.307	0.128	0.090	4008.0	2905.0	1.0
Overall_mean[12]	511.985	7.413	496.975	524.942	0.128	0.091	3351.0	3018.0	1.0
Overall_mean[13]	493.300	11.085	473.001	514.334	0.181	0.128	3741.0	2519.0	1.0
Overall_mean[14]	503.008	6.023	491.690	514.593	0.099	0.070	3700.0	3425.0	1.0
Overall_mean[15]	495.492	7.078	481.897	508.810	0.120	0.085	3512.0	3296.0	1.0
Overall_mean[16]	515.685	7.168	501.659	528.930	0.111	0.078	4167.0	3459.0	1.0
Overall_mean[17]	492.426	6.316	481.276	505.205	0.106	0.075	3553.0	3062.0	1.0
Overall_mean[18]	477.580	8.455	461.670	493.446	0.140	0.099	3650.0	2687.0	1.0
Overall_mean[19]	508.198	7.648	494.250	523.141	0.129	0.092	3486.0	3499.0	1.0
Overall_mean[20]	354.496	13.938	329.092	382.214	0.204	0.144	4657.0	3581.0	1.0
Overall_mean[21]	520.310	7.693	506.408	535.676	0.119	0.084	4151.0	3656.0	1.0
Overall_mean[22]	501.337	8.565	486.377	518.502	0.155	0.110	3048.0	2868.0	1.0
Overall_mean[23]	489.883	8.175	475.078	505.824	0.143	0.101	3297.0	2963.0	1.0
Overall_mean[24]	501.974	9.961	483.941	521.693	0.181	0.128	3056.0	3390.0	1.0
Overall_mean[25]	403.136	11.231	381.858	424.031	0.195	0.138	3296.0	3159.0	1.0
Overall_mean[26]	381.531	17.842	349.633	417.236	0.292	0.207	3729.0	3202.0	1.0
Overall_mean[27]	472.908	10.153	453.432	491.716	0.168	0.119	3642.0	2786.0	1.0
Overall_mean[28]	477.533	11.546	455.395	498.945	0.195	0.138	3497.0	2746.0	1.0
Overall_mean[29]	489.548	12.874	465.201	513.021	0.257	0.182	2538.0	2853.0	1.0
Overall_mean[30]	431.785	7.851	418.175	447.797	0.119	0.085	4328.0	2916.0	1.0
Overall_mean[31]	449.827	7.590	436.179	465.141	0.121	0.086	3925.0	2895.0	1.0
Overall_mean[32]	512.782	6.730	500.580	526.278	0.111	0.079	3672.0	3705.0	1.0
Overall_mean[33]	521.225	11.537	499.554	542.793	0.208	0.147	3089.0	2917.0	1.0
Overall_mean[34]	540.586	9.612	522.033	558.469	0.153	0.108	3957.0	3531.0	1.0
Overall_mean[35]	464.240	5.331	454.691	474.984	0.076	0.054	4970.0	3029.0	1.0
Overall_mean[36]	450.390	8.991	433.623	467.127	0.144	0.103	3836.0	2568.0	1.0
Overall_mean[37]	465.725	19.685	429.067	502.430	0.393	0.278	2541.0	2519.0	1.0
Overall_mean[38]	448.849	7.163	435.301	462.651	0.104	0.073	4774.0	3376.0	1.0
Overall_mean[39]	447.518	18.718	411.696	481.451	0.361	0.256	2693.0	2576.0	1.0
Overall_mean[40]	441.304	7.067	428.375	455.446	0.108	0.076	4310.0	3244.0	1.0
Overall_mean[41]	480.386	7.224	466.706	493.810	0.112	0.079	4156.0	2848.0	1.0
Overall_mean[42]	494.988	5.671	484.411	505.879	0.086	0.061	4349.0	3137.0	1.0
Overall_mean[43]	489.515	7.013	476.526	502.930	0.120	0.085	3441.0	3267.0	1.0
Overall_mean[44]	511.259	7.157	498.045	525.311	0.117	0.083	3752.0	3045.0	1.0
Overall_mean[45]	504.267	7.923	488.608	518.724	0.132	0.094	3599.0	3164.0	1.0
Overall_mean[46]	386.966	12.137	363.771	409.502	0.185	0.131	4294.0	3153.0	1.0
Overall_mean[47]	395.127	11.023	375.705	416.690	0.173	0.122	4087.0	3388.0	1.0
Overall_mean[48]	398.730	13.489	372.901	424.018	0.235	0.166	3346.0	2641.0	1.0
Overall_mean[49]	485.227	7.337	471.452	498.994	0.130	0.092	3187.0	2799.0	1.0
Overall_mean[50]	522.528	9.670	503.077	539.189	0.181	0.128	2859.0	2608.0	1.0
Overall_mean[51]	438.286	8.154	423.132	453.733	0.145	0.102	3178.0	2889.0	1.0

since it’s less clear what the spline coefficient terms mean.

az.summary(spline_idata)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
Intercept	225.689	205.515	-133.222	630.021	6.036	4.292	1163.0	1653.0	1.0
bs(Income, degree=3, knots=knots_income)[0]	-355.159	1613.343	-3460.788	2642.535	56.357	39.865	819.0	1504.0	1.0
bs(Income, degree=3, knots=knots_income)[1]	194.905	122.719	-42.478	423.960	2.343	1.728	2748.0	2719.0	1.0
bs(Income, degree=3, knots=knots_income)[2]	48.748	123.638	-182.012	272.863	4.065	2.875	923.0	1431.0	1.0
bs(Income, degree=3, knots=knots_income)[3]	68.533	98.242	-118.858	250.496	3.369	2.383	853.0	1416.0	1.0
...	...	...	...	...	...	...	...	...	...
Overall_mean[47]	391.621	17.972	358.809	425.230	0.307	0.218	3420.0	3533.0	1.0
Overall_mean[48]	433.747	21.807	392.182	474.390	0.346	0.245	3955.0	3931.0	1.0
Overall_mean[49]	527.171	14.085	498.983	551.381	0.216	0.153	4288.0	3655.0	1.0
Overall_mean[50]	512.424	17.480	480.801	546.508	0.256	0.181	4656.0	3790.0	1.0
Overall_mean[51]	412.105	14.369	385.225	439.530	0.221	0.156	4238.0	3730.0	1.0

78 rows × 9 columns

We can also check that the model seems to recover the observed data well.

fig, axs = plt.subplots(2, 1, figsize=(8, 12))
axs = axs.flatten()
az.plot_ppc(base_idata, ax=axs[0])
az.plot_ppc(spline_idata, ax=axs[1])

array([<Axes: xlabel='Overall'>], dtype=object)

Next we’ll dig into the spline basis features and decompose the predicted outcome.

The Component Contributions

This will also serve as a kind of posterior predictive check for each country as we can show the degree which posterior draws from each component sum to achieve a plausible mirror of the observed data.

Bincome = spline_model.response_component.design.common['bs(Income, degree=3, knots=knots_income)']

income_coefs = az.extract(spline_idata['posterior']['bs(Income, degree=3, knots=knots_income)'])['bs(Income, degree=3, knots=knots_income)']

Bedu = spline_model.response_component.design.common['bs(Edu, degree=3, knots=knots_edu)']

edu_coefs = az.extract(spline_idata['posterior']['bs(Edu, degree=3, knots=knots_edu)'])['bs(Edu, degree=3, knots=knots_edu)']


Bhealth = spline_model.response_component.design.common['bs(Health, degree=3, knots=knots_health)']

health_coefs = az.extract(spline_idata['posterior']['bs(Health, degree=3, knots=knots_health)'])['bs(Health, degree=3, knots=knots_health)']

income = np.dot(Bincome, income_coefs).T 
edu = np.dot(Bedu, edu_coefs).T
health = np.dot(Bhealth, health_coefs).T

intercept = az.extract(spline_idata['posterior']['Intercept'])['Intercept'].values

fig, ax = plt.subplots(figsize=(9, 7))
for i in range(100):
    if i == 1:
        ax.plot(income[i], label='Income Component', color='red')
        ax.plot(edu[i], label='Edu Component', color='blue')
        ax.plot(health[i], label='Health Component', color='darkgreen')
        ax.plot(intercept[i] + income[i] + edu[i] + health[i], label='Combined Components', color='purple')
    else: 
        ax.plot(income[i], alpha=0.1, color='red')
        ax.plot(edu[i], alpha=0.1, color='blue')
        ax.plot(health[i], alpha=0.1, color='darkgreen')
        ax.plot(intercept[i] + income[i] + edu[i] + health[i], color='purple', alpha=0.3)

ax.scatter(range(len(spline_idata['observed_data']['Overall'])), spline_idata['observed_data']['Overall'], label='Observed', color='grey', s=56, ec='black')
ax.set_title("Additive Spline Components", fontsize=20)
ax.legend();
ax.set_xticklabels(pisa_df.dropna(axis=0).reset_index()['Country'])

[Text(-10.0, 0, 'Argentina'),
 Text(0.0, 0, 'Australia'),
 Text(10.0, 0, 'Austria'),
 Text(20.0, 0, 'Belgium'),
 Text(30.0, 0, 'Brazil'),
 Text(40.0, 0, 'Bulgaria'),
 Text(50.0, 0, 'Canada'),
 Text(60.0, 0, 'Chile')]

We can see here how the combined components borrow the structure of the outcome variable primarily from the income variable. The health measures have closter to zero additive contribution while the uncertainty in the educational component varies wildly. But the blue educational component here is used primarily as a scaling contribution which adds to the level of the outcome variable rather than distorting the shape.

Conclusion

We’ve seen the application of splines as univariate smoothers to approximate wiggly curves of arbitrary shape. We’ve also tried gaussian process approximations of the same univariate functions. The suggested flexibility of both methods is a strength, but we need to be careful where splines have a tendency to over-fit to individual curves. As such we have tried to show that we incorporate spline basis modelling in a hierarchical bayesian model and recover very compelling posterior predictive checks and additionally derive predictiions from the mixed variant of the hierarhical model which helps us understand the implications of the data for generic forecasts of insurance loss curves. Finally we showed how splines can be used additively to model non-linear functions of multiple covariates. These are a powerful tool to interrogate complex non-linear relationships.